{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 旅行商问题\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览\n",
    "\n",
    "旅行商问题（travelling salesman problem, TSP）是组合优化中最经典的 NP–困难的问题之一，它指的是以下这个问题：\"已知一系列的城市和它们之间的距离，这个旅行商想造访所有城市各一次，并最后返回出发地，求他的最短路线规划。\"\n",
    "\n",
    "这个问题也可以用图论的语言来描述。已知一个有权重的完全图 $G = (V,E)$。它的每个顶点 $i \\in V$ 都对应一个城市 $i$，并且每一条边 $(i,j) \\in E$ 的权重 $w_{i,j}$ 对应城市 $i$ 和城市 $j$ 的距离。需要注意的是，$G$ 是个无向图，所以权重是对称的，即 $w_{i,j}=  w_{j,i}$。根据以上定义，旅行商问题可以转化为找这个图中最短的哈密顿回路（Hamiltonian cycle）的问题。哈密顿回路为一个通过且仅通过每一个顶点一次的回路。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用量子神经网络解旅行商问题\n",
    "\n",
    "使用量子神经网络方法解决旅行商问题，需要首先将该问题编码为量子形式。\n",
    "具体的包括以下两部分：\n",
    "\n",
    "1. 旅行商经过城市的顺序将编码进量子态 —— ${\\rm qubit}_{i,t} = |1\\rangle$ 对应于在时间 $t$ 经过城市$i$。\n",
    "    1. 以两城市$\\{A,B\\}$举例。先经过$A$再经过$B$ 将由$|1001\\rangle$表示，对应于旅行商在时间 $1$ 经过城市$A$，在时间 $2$ 经过城市$B$。\n",
    "    2. 类似的 $|0110\\rangle$对应于先经过$B$再经过$A$.\n",
    "    3. 注意：$|0101\\rangle$意味着在时间 $2$ 同时经过城市 $A$、$B$，而这是不可能的。为避免此类量子态，我们会通过引入代价函数的方式 (具体见下一节)。\n",
    "\n",
    "2. 总距离被编码进损失函数： \n",
    "\n",
    "$$\n",
    "L(\\psi(\\vec{\\theta})) = \\langle\\psi(\\vec{\\theta})|H_C|\\psi(\\vec{\\theta})\\rangle \\, ,\n",
    "\\tag{1}\n",
    "$$\n",
    "其中 $|\\psi(\\vec{\\theta})\\rangle$ 对应于参数化量子电路的输出。\n",
    "\n",
    "在下一节中将详细介绍如何编码旅行商问题为对应量子问题。通过优化损失函数，我们将得到对应最优量子态。再通过解码，将得到最终的路线规划"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 编码旅行商问题\n",
    "\n",
    "为了将旅行商问题转化成一个参数化量子电路（parameterized quantum circuits, PQC）可解的问题，我们首先需要编码旅行商问题的哈密顿量。\n",
    "\n",
    "我们先将此问题转化为一个整数规划问题：假设图形 $G$ 的顶点数量为 $|V| = n$ 个，那么对于每个顶点 $i \\in V$，我们定义 $n$ 个二进制变量 $x_{i,t}$，$t \\in [0,n-1]$：\n",
    "\n",
    "$$\n",
    "x_{i, t}=\n",
    "\\begin{cases}\n",
    "1, & \\text{如果在最后的哈密顿回路中，顶点 } i \\text { 的顺序为 $t$，即在时间 $t$ 的时候被旅行商访问}\\\\\n",
    "0, & \\text{其他情况}\n",
    "\\end{cases}.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "因为 $G$ 有 $n$ 个顶点，所以我们共有 $n^2$ 个变量 $x_{i,t}$，所有这些变量的取值我们用 $x=x_{1,1}x_{1,2}\\dots x_{n,n}$ 来表示。现在我们假设 $x$ 对应了一条哈密顿回路，那么对于图中的每一条边 $(i,j,w_{i,j}) \\in E$，条件 $x_{i,t} = x_{j,t+1} = 1$（也可以等价地写成 $x_{i,t}\\cdot x_{j,t+1} = 1$）成立当且仅当该哈密顿回路中的第 $t$ 个顶点是顶点 $i$ 并且第 $t+1$ 个顶点是顶点 $j$；否则，$x_{i,t}\\cdot x_{j,t+1} = 0$。所以我们可以用下式计算哈密顿回路的长度\n",
    "\n",
    "$$\n",
    "D(x) = \\sum_{i,j} w_{i,j} \\sum_{t} x_{i,t} x_{j,t+1}.\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "根据哈密顿回路的定义，$x$ 如果对应一条哈密顿回路需要满足如下的限制：\n",
    "\n",
    "$$\n",
    "\\sum_t x_{i,t} = 1 \\quad  \\forall i \\in [0,n-1] \\quad \\text{ 和 } \\quad \\sum_i x_{i,t} = 1 \\quad  \\forall t \\in [0,n-1],\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "其中第一个式子用来保证找到的 $x$ 所代表的路线中每个顶点仅出现一次，第二个式子保证在每个时间只有一个顶点可以出现。这两个式子共同保证了参数化量子电路找到的 $x$ 是个哈密顿回路。所以，我们可以定义在此限制下的代价函数：\n",
    "\n",
    "$$\n",
    "C(x) = D(x)+ A\\left( \\sum_{i} \\left(1-\\sum_t  x_{i,t}\\right)^2 +  \\sum_{t} \\left(1-\\sum_i  x_{i,t}\\right)^2  \\right).\n",
    "\\tag{5}\n",
    "$$\n",
    "\n",
    "其中 $A$ 是惩罚参数，它保证了上述的限制被遵守。因为我们想要在找哈密顿回路的长度 $D(x)$ 的最小值的同时保证 $x$ 确实表示一个哈密顿回路，所以我们需要设置一个大一点的 $A$，最起码大过图 $G$ 中边的最大的权重，从而保证不遵守限制的路线不会成为最终的路线。\n",
    "\n",
    "我们现在需要将代价函数 $C(x)$ 转化为一个哈密顿量从而完成旅行商问题的编码。每一个二进制变量可以取0和1两个值，分别对应量子态 $|0\\rangle$ 和 $|1\\rangle$。**每个二进制变量都对应一个量子比特，所以我们需要 $n^2$ 个量子比特来解决旅行商问题。** 和最大割问题相似，因为泡利 $Z$ 的两个本征态和我们需要的量子态 $|0\\rangle$、$|1\\rangle$ 一样，这两个本征态所对应的本征值分别是 1 和 -1，所以我们考虑利用泡利 $Z$ 矩阵将代价函数编码为哈密顿量。\n",
    "\n",
    "我们现在将二进制变量映射到泡利 $Z$ 矩阵上，从而使 $C(x)$ 转化成哈密顿矩阵：\n",
    "\n",
    "$$\n",
    "x_{i,t} \\mapsto \\frac{I-Z_{i,t}}{2}, \\tag{6}\n",
    "$$\n",
    "\n",
    "这里 $Z_{i,t} = I \\otimes I \\otimes \\ldots \\otimes Z \\otimes \\ldots \\otimes I$，也就是说 $Z$ 作用在位置在 $(i,t)$ 的量子比特上。通过这个映射，如果一个编号为 $(i,t)$ 的量子比特的量子态为 $|1\\rangle$，那么对应的二进制变量的取值为 $x_{i,t} |1\\rangle = \\frac{I-Z_{i,t}}{2} |1\\rangle = 1 |1\\rangle$，也就是说顶点 $i$ 在最短哈密顿回路中的位置是 $t$。同样地，对于量子态为 $|0\\rangle$的量子比特 $(i,t)$，它所对应的二进制变量的取值为 $x_{i,t} |0\\rangle = \\frac{I-Z_{i,t}}{2} |0\\rangle = 0 |0\\rangle$。\n",
    "\n",
    "我们用上述映射将 $C(x)$ 转化成量子比特数为 $n^2$ 的系统的哈密顿矩阵 $H_C$，从而实现了旅行商问题的量子化。这个哈密顿矩阵 $H_C$ 的基态即为旅行商问题的最优解。在接下来的章节中，我们将展示怎么用参数化量子电路找到这个矩阵的基态，即对应最小本征值的本征态。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paddle Quantum 实现\n",
    "\n",
    "要在量桨上实现用参数化量子电路解决旅行商问题，首先要做的便是加载需要用到的包。其中 `networkx` 包可以帮助我们方便地处理图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:15.901429Z",
     "start_time": "2021-05-17T08:00:12.708945Z"
    }
   },
   "outputs": [],
   "source": [
    "# 加载量桨、飞桨的相关模块\n",
    "import paddle\n",
    "import paddle_quantum\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum import Hamiltonian\n",
    "\n",
    "# 旅行商问题相关函数\n",
    "from paddle_quantum.QAOA.tsp import tsp_hamiltonian  # 构造旅行商问题哈密顿量的函数\n",
    "from paddle_quantum.QAOA.tsp import solve_tsp_brute_force  # 暴力求解旅行商问题\n",
    "\n",
    "# 用于生成图\n",
    "import networkx as nx\n",
    "\n",
    "# 加载额外需要用到的包\n",
    "from numpy import pi as PI\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 生成该旅行商问题中的图 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们生成该旅行商问题中的图 $G$。为了运算方便，图中的顶点从0开始计数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.212260Z",
     "start_time": "2021-05-17T08:00:15.918792Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0UElEQVR4nO2deXhV5bX/PysTSUAGEUUUiQOg4IQiggKiOEeL81RpxWptq7bVn1fN5Wqt1xovtVJrr1rFqx2ciqW2GgccKU6IOBQZRNEggzMyhJzM6/fH2oecDEACJ2effbI+z3MenuTss/ciyfmetde73u8SVcVxHMdJDVlhB+A4jtOZcNF1HMdJIS66juM4KcRF13EcJ4W46DqO46QQF13HcZwU4qLrOI6TQlx0HcdxUoiLruM4Tgpx0XUcx0khLrqO4zgpxEXXcRwnhbjoOo7jpBAXXcdxnBTious4jpNCXHQdx3FSiIuu4zhOCnHRdRzHSSEuuo7jOCnERddxHCeF5IQdgOM4qaWopCwL2BMYBBQAeUANEAOWAEvLS4sbwoswsxGfBuw4mU0gsuOBYmAMsA/QANQBEjw0eORgd8CLgNlAGfCCi3DycNF1nAylqKSsFzAJuArYDuiKCWxbUWADsA74DXB/eWnxt8mOs7Phous4GUZRSVkh8D/ARVhGW5iE01ZiGfA04Jry0uLKJJyzU+Ki6zgZRFFJ2RjgEaAXVq9NNjHgW+Ds8tLiVzrg/BmPi67jZABFJWVdgKnABXSM2DYnBjwAXFFeWlydgutlDC66jhNxikrKugEzgQNJjeDGiQHvAMeVlxZXpPC6kcZF13EiTCC4rwCDgfwQQqgCPgBGu/C2Dd8c4TgRJSgpzCQ8wSW47mDg2SAeZwu46DpOdJmKlRTCEtw4+cAw4LaQ44gEXl5wnAgSdCk8S2pruFsiBhzrXQ2bx0XXcSJG0If7IdAv7FhaYRUw0Pt4N42XFxwnekzB+nDTkV7ALWEHkc54pus4ESLY2ruK8Ou4m6MK6OdbhlvHM13HiRaTsK296UwDtknDaQXPdB0nIgRuYSuAncOOpQ2sAvq7O1lLPNN1nOgwHnMLiwLdgaPCDiIdcdF1nOhQjNkzRoFC4MSwg0hHXHQdJzqMoX1+uGGSBYwNO4h0xGu6jhMBgnruBpLUtdAlJ4uSE/bhpP13pluXHN5fuZabnlrEu8vXJOP0cWJA1/LSYheZBDzTdZxosCdQn6yTXX/SEC44rIivK6qZufBzDtqtF3++cAS9CnOTdQmwLoY9k3nCTMBF13GiwSBsptk207trHmce3J/6BuW70+bw00fe5fH3VrJdfi7fH1WUjEvEqcPidhJw0XWcaFBAkuq5g3bajrycLFatifHNhhoA5q9YC8CQnbsn4xJxhPTyhkgLXHQdJxrkkSTR3aFbHgAbahoT58oaq1z02S6p7owCuN1jM1x0HSca1GDTebeZryssu+2al7Pxe127ZAPw1fqkTt5RwEf5NMNF13GiQYwkie6HX66npq6Bfj0LNma9++/aE4BFn69LxiXiKBa3k0DOlg9xHCcNWKKqOSLbXmH4uqKGx95ewXkjduPBH4xkyRfrKd5vZyqq6/jj68uSEOpGcoAlyTxhJuCi6zhpioj0Aw4FDkVkRP8rH+squckpkf7yiQXU1TdQvN/OFPXeiXeWr+FXTy1kdbCwliSygKXJPGEm4JsjHCcNEJGuwMHERdYeuyYe0/eC39Kl714hRLfVzCsvLR4edhDphme6jpNiRCQb2JumArsvkN3s0PXAm8Ac4M2cHn2LgYuIxlbgBuBfYQeRjrjoOk4HIyI701Rgh9PSLaweeBcT2PhjsaputEYsKimrBM4FunV81NtMJfBU2EGkIy66jpNERKSQlmWC/q0c+imNWewcYJ6qbmmu2AtY9hsF0V0HvBh2EOmI13QdZysRkSxgHxrFdQSwH62XCebSKLBvqupnW3PNopKyK4H/xqwT05VK4L/KS4unhh1IOuKZruO0ERHpS9MM9hA2XSZIzGIXq2qyzGruB36VpHN1FFnAA2EHka54puukC9lAT6yZPvTx3QllghE0iuxurRy6nKZ12LdVdUNHxlZUUnYH8APS09cgBkwrLy3+adiBpCue6TphIFg71KHA4dhYl72xFe8a4EfAIyRpB9YWg7EyQfNugtbKBBU0LRPM2doywTZyDXAa6Sm6q4Frww4infFM10kVAkwAfgiMwsy4a7BFoebb0SuBEuB3HRKIyE60LBM0t9dqAObTspsgaZ6220JRSdloYCbpJbwx4Jjy0uJXww4knXHRdVLFD4HbaPuMr0+APbb1okGZ4CCaimxrZYIVNBXYeR1dJthWikrK7sRGnaeD8MaA+8tLiy8NO5B0x8sLTiroS/sEF6z80Av4tq0vCMoEg2kqsPuz+TLBm1iZYFU7YksXrgAOwD5UkjLGZyupAt4Brgwxhsjgma6TCo4BpgM9Er/5zTffMGfOHGbPns1ee+3FBRdcQHb2Rn1cC5yF3UK3SrMywQisTNCj2WENwPs0zWIXpUuZYFspKinrpqqvUF+3n+TkhuEaWAV8AIwuLy2uCOH6kcMzXScVfEErNqKXX345a9euZfTo0ZSVlVFdXc2Pf/xjAietQqz2OxNARApoWSYY0Mq14mWCeMvWPFXNWDEoLy2uyC7s8XKf0687IG+nPchKkiFOG4lhGe5xLrhtp9NmusF01T2xGU4FmDN/DfaHtARYWl5a3LDpMzjtIBe7pc9L/GZDQwNZWabFd9xxB3PnzuX+++/fmO2uWLFiSf/+/V+gsUzQPEnYQMtugiiWCbYaEbkAuJ/s3Lp+F97xZG7vXY8jNTXeGNYzfGV5abEblbeDTiO6gciOB4qBMdhOogZseJ4EDw0eOVhmtgiYDZQBL7gIbxPzMVOXJnz55ZfcfvvtzJ49m4svvpiJEydufG79+vV0776xqSBeJkjcdLAwU8oEW4OIHI5ttc0DLlHVe4KuhkexenhHiG8Mq7OfXV5a/EoHnD/jyXjRLSop6wVMAq7Cdg91pX0uTYplVOuA32ArtG1e3HE28jvgMpr97L/66ituuOEGdt99d55//nluv/12Bg8eDEBVVVX98ccff+usWbOeJsPLBO1FRIqwD6A+wO9U9Wfx54pKygqBW4CLsQ+rZGwZrsQSkXuBa8tLi0PfwBJVMlZ0gz+8/8Gs8JL9hzcNuMb/8DZP0E0wCDj06quvPu+6664b361bt+adBBu59dZbWb9+Pb/4xS/iZYf1wCXAw6mJOBqISDfgNWwDx0ygWFVbjGcPEo4LsISjO/YeaPNimzY0IFlZFVjCcSvwgCcc205Gim5RSdkYbEeT32KlEBHpQ8tNBz0BdtttNxYtWkRh4aY/+6677jpqa2u5+eab46KrwN3ATzo49MgQfJDNwDaaLAEOVdU1m3tNQmntBGAsMITNlNZUNavmi4/zqj+dn53bZ8DEgt2HPeSlteSRUaJbVFLWBZhK6hrGY5ixxxWdbTFBRPJp7CaI+xPs3sqhq4A5WVlZc6qrq6/PycnZqLq1tbV89NFHzJ49m8WLF/Pee+9x/fXXc8QRRyS+fiEwtOP+J9FCRG7GduutwQS33TPIAhHeg8ZF5C7Y1N6Ni8jLbjnpUeBM4FJVvTM50TuQQaJbVFLWDbvVOpDU7tDJ+LaZILsaSNMs9gBa7yZ4i6YWhisSnn8JGJf4gkceeYQZM2YwZswYxo0bx3777df88nXYQlFm/KFuAyJyPvBnzMnseFV9vgOvdSFwH/CEqn6no67TGckI0Q0E9xVsN1IYO3MyqkE8oUwwIuHfns0OU2ABTTcdLGyttpjAZOAXWAtZW6kDemN1xU6LiIwEXsay0stU9X87+Hq7YD3PG4Deqtqp7uQ6ksiLblBSeAkYRvhbId8GjopSqSEoEwyjaRa7yTIBjRsP3lLV9e283Hjgb7TcNbYp4qYzw+jEma6I9Mf6kXcC7lLVlNS4ReQ9rD96vKr6FIgkkQk70qZiJYUwBZfg+sMwj4G0NP0IygR70bJM0DzzrKRpmWBOszLB1jIHy9TaQm3wOIfOLbhdgX9ggvsi8LPNvyKpPIOJ7vH46J2kEelMN+hSeJb0cFmKEwOOTYeuBhHZgaYCu6kywUKalgkWbKFMsC38ETidluY3MUxkC4APgaex1rzFHRRH2hN8SP4V+3ktBUao6uoUXv9ITGzfV9UWxXZn64is6AZ9uB8C/cKOpRVWAQNT2ccblAkOpKnItmaN+BlNBXZrygTbQhZmUn4rlvV+AbyOlYjmAO9h27E7PSLyS+B6rJ49UlUXpfj6eZgpeVegf5Ludjo9US4vTMH6cNORXtiOoA4ZWSLmCJPYTTACE9zWygTzaCqyKzTcT9oG4E7gHmwVPpqf+h2MiJyNCW4DcHaqBRdAVWtE5AXgO8BxWDeDs41EMtMNdtqsIvw67uaoAvolYwdPUCZInNU1gpYfOKkuEzgdhIgMxzw/8oErVPW3IcbyY+xD8jFVPTOsODKJqGa6k7AMIJ1pwDZptGsMtYh0oWWZYM9WDv2clmWCTt1WlQkErVr/wAR3GnB7uBHxTPDvMSKS4x/i207kMt1gN80KYOewY2kDq4D+m9pCGZQJmncTHEjLMkGMZpsOgOUhlwmcJBN4Bv8LGB78e4yqhl7fFpHFWA/8GFUNfYE46kQx0x2PuYVFge7YpNvnAUSkNy3LBNs3e01rZYL3PcPIbIIP4P/DBLccOD0dBDfgGUx0j8c2ITnbQBRFt5j2zdoKDVUtrPp0/mSRkyZhArtXK4fFywRxn9i5XibolEzGepIrgJNV9euQ40nkGaw/+Hjgv0KOJfJEUXTH0D4/3NAQkaysvIJxCd+K0bKbwMsEnRwROR34b+wu51xVfT/kkJozC1sYPlhEdlTVL7F+76ODx16Yz+5f8W6ULRIp0Q3quUOSdb4LDy/izIP7M2in7cjOEn77/BJ++8KHyTo9AHk7FtUjchmq8TJBbVIv4EQaERkG/Cn48hpVfTLMeFpDVWPZ2dmzhg0bdlxpaemdmDvZYEyIt8OSoEOBa4DDgu87myBSoout4idtPMu+u/RgbayWz9bG2LVXMjzOWyLZuVUDrnni+fLS4o865AJOZBGRvsA/MXPxP2IbRtKJPlh/7mmxWGxcVVUVBQUFp9JohJ44864b1jt+LjY7zdkEURPdQZjrVFK48q/vAXDP+Qd3mOhi8Q4CXHSdjQQ7CP8O7IpNgbgkjcpMpwK/xP5ua4Dt8vLyyMvLg81PnuiGjQhy0d0MURPdAiJSz01ASC9vCCdkgk6Fe4GRwKfAaWlknTgK+AuN463aO9P9YOzvPZbMoDKJNs9LShPyiKbotvcP18lsrgbOx7xqv6OqX4QcTyJ/YCuSBFVl7dq1YGJ7xBYO79RETXRriN7qqGKjUBwHEfkOUBp8eb6qvhdmPM0QbIGsXYlNfX09TzzxBBMmTABbWPNJE5shaqIbI5qi67daDiKyP/AQJmqTVfXxcCNqgQLLN3fAypUrW3yvoaGBvffem7lz5/Luu+9m4aK7WaJW011CEmM+e3h/DinqxdBdbJDBsUN2YtdeBcxc+AUzFybtji8Hi9vpxIjIjlinQldMeEs3/4rQ+Cfmjpfd2pMTJ05kypQp7L///jz33HO8/PLLzJkzh4qKCo4++mjy8/PBdlnuDnySurCjQ9REdylJzM4PKerFGQf33/j1kH49GNKvByu+jSVNdFU1W0SWJuVkTiQJTIz+BgzAdh5elEadCs15AvgBtoW9BcOGDaO4uJjdd9+drKwsRo4cyfXXX89RRx1FVtbGt2YDtnvtrtSEHC2iaHgzDxv9HQmqP/uQz/94xVLMwCT++CSN33ROEgk6Fe7DnPFWAoeo6mfhRrVZ8oC1bMI2dfbs2Zx44ol89dVX8ax2IytXrmSXXXaJf/kS5jviNCNqmS6Yz+gwItDFoA0NWr1yUS22qWNP7I0HsFJEEkV4kYtwxnIF9nuPYZ0K6Sy4YIvVrwNHtvbkmDFjmDBhAvn5+axZs4ZHH32Uv/3tbyxdupQDDzyQQw89lOLiYoYOHToKE/B0Me1JG6KY6R4DzMAasdOdiobaqtOX/+aMb4CxwWMMNlI8ka+xD5N/Yfvc/62qSdt554SDiJyI3a5nAWep6vSQQ2orP8Z2x7W6Y2jNmjW8/vrr3HzzzXTp0oWTTz6ZCRMmUFtbyx133MG6det44IEH1gGnYBmvk0AUM90XgPVEQ3TXZeXmP6+qDZjRzdRg2OA+NIrwEZg38KnBA2CdiLxCYyY8L41s/pw2ICJDgEcwwb0hQoIL5ir2m009WVFRwYMPPsjVV1/NySef3OS58847j8svvxxswfBkXHRbELlMF6CopOxKzJWpw/buJoFK4L/KS4s3OzkiqPntQaMAj8VWfhOJYbd8cRGeo6opG3rptI9gvNIc7Pf6V+CcCJaPVgC7tPbEqlWrGDduHEuWNG3KWbZsGVOnTuXggw9m4sSJYN0LrQ1H7dREVXQzekaaiPTHyhDxbHifZofUAnNpLEe85h686UEwQXcm9gE6Dxgb0Q/IO4FL2ES30NixYxk+fDijRo1i3rx5lJWV8e2333L++eczefJktttuO7C/0z7YwpwTEEnRBSgqKbsDa21JR1+DGDCtvLQ4KdOAgx7P0TRmwwfQdCGxAXiHxkz4lTQzwe4UBHctdwM/xEbdH6KqLXcTRIPjgUfZROvYxx9/zFtvvcXjjz9O9+7dmTBhAieccELzw2LYovcHHRtqtIiy6BYCHwL9wo6lFVYCg8pLizskwxGRnsDhNGbCw2lZn19AQpuaqq7qiFicRkTkcuB32F3OEar6ZsghbQsFwLe00zekoaEBEUFEFPg3JrrRFJkOIrKiC1BUUjYau5VLp2w3BhxTXlr8aqouKCJdMcequAiPpGXp5SOa9gqXR7DOmLaIyLHA09jt+Hmq+nDIISWDVzFT8s2SILTxb1Vjd19HYrVtJ4FIiy5AUUnZndio83QQ3hhwf3lp8aVhBhHsgBpOYznicFp2e6ygqQgvdhHeOkRkMCYuPYCbVPW6kENKFj/Ftiu3ZcF6PZYVLwWmY25lfnfVCpkgul2AF7FdamEurFUBbwNHlZcWp5WrmIjkYKPdE3uFm08h/grrFZ6FifB87xXeMiKyPfAGNjXh78AZQYtgJrADZoDT2vuqura2VisrK/Pnzp37+dFHH30F8BzwTUojjCCRF12AopKybtho6MGEI7xV2GLB6PLS4ooQrt8ugl7hITTtFe7b7LC12M80sVfY57slICK5WElhPPAecLiqbgg3qqRzMma4XoiVDAqAhcD0559//o1jjz32BVWtArZXVXfTawMZIbqwUXhnYhldykoNqhoTkXeA46IguK0RrLrvSWOf8FigqNlhlbTsFU71m2w37Pe7HJhPEkc3bQ0i8nvgUuALYISqfhpmPB1IF+AkbGFtHgktYCLyFjYt4gRVfSac8KJFxogubCw13Ibtde9w4W2orabm849m5/cfeky6lRS2FRHZjaa9wns3O6QGc8yKi/Brqrq+g8LJAUqCRy3WLpeLtcmVYSWRuaTQLF5Efoz1stYA41T19VRdO50QkZuAycDtqvrzkMOJBBklunGCroZHgV50jPjGGmqrq7589Lpe1SsW1gCjVPXtDrhO2iAiO9G0V3h/WvYKv03TXuFk1fcmAb+n9QWdGqy8k49lv3cBf6YDjVZE5Cjsriob+J6q/rmjrpXuiMhobC3gA1Vt/sHstEJGii5s7OO9BZtO2kBytgxXYi1B9wLXLrvlpFsxc5CPgIM6MNNLO0SkFy17hZsbX79P017hrXXYit/CtoVK4B7M3WtryMF+x62KtogMxDoVegH/o6rXbuV1MoJgkfZrrHNjD1V14/ItkLGiGyfYMnwBcBW2u6aQ9hmhN2Bv5HWY89ID8a29wRjtOVjW9yAwsbO2XYlIN1r2CjdvrP+Qpm1qy9rw8+qC/ezz2hFODPsQWNjG4wW4FitfdMV+3z8H/o+Exv5gU8ob2ILtE8Cp3uEBIjIdOAP4iaq6cfkWyHjRjVNUUpaFrTKfgInCEExQ67A3nWBvMKUx21mIicNTwIvlpcUtWoFEZG9scaEQmKSqD3T0/yUKBL3Ch9C0V7hrs8OW01SEP2hFhMcAT9JsO+qGDRuYOXMm++67LwMHDmx++Q3A5cD9bQz3OuCaZvFtwKY9XAzUBBndk8BxWAZ/WGe6s9kcIvIDYBrwT1WdEHY86U6nEd3mBCK8BzAIq/t2wRZiYthMs6XlpcVt+uGIyAXYG7wSOFhVF3dEzFEmEK1hNO0V7tXssK9oFOBZwPuqWgJcjy2cbWT+/PncdNNNfPzxxwBMnTqV0aNHx59ei2Vez7chtFOBv9B6+akS++A9QUQmY9nv15inQnkbzt0pEJFdsQ/QDVjrmNuQboZOK7rJJGi5+jPwXWy/+aFB76KzCYJe4aE07RXeqdlhaxYsWKBDhgxpLs7EYjEKCmyNdNq0acyaNYs//3njelYt1ti/Jee1A4DX2Hy9v2bDhg2xESNG9Fi4cGEtMF5VZ2/hvJ0OEZkP7AscparuobsZojaCPS0JbonjC2r7sxkDaMdQ1QZVna+q/6uqZ2NG7oOx2/k/A8uys7N77rHHHi0EF6CgoICaGkuoqqur6datG7W1G/dulLNlwd0J20G1pe6WvIKCgh5z5szhpptuutMFd5M8HfzbwmrMaYaq+iNJD2wrcg1WFz4t7Hii/pgxY8bJ1dXVldoKNTU1qqq6evVqHTZsmD788MPxp+pV9Y4tnLuLqr6rqjWJ51y9enVrl0qkUlWvUVUJ+2eTbg9sCKVio6ZCjyedH57pJhG1Xt3/CL68T0SKQgwn8px66qmD8vLyWh1Ampuby7PPPsuRRx7JpEmTOOeccwCoqKjQq666qr+IfCfwRWiNX2O1/I114iuuuIILLriAs846i7lz524qpAJs0e0h2ml52Al4Favp7icirU6ccALCVv1Me2BdEP/EPvVfA3LDjinCj5e0FVauXKlTpkzRM844Q2fOnNnkuaqqKu3du3e8CyXu6fp74Eygr6r20GYZ7s0336xHH320fvvtt3rttdfqMccco7FYrLVLx6lU1dmqmpcGP6O0eST83V8Ydizp/Ag9gEx8YNN+lwd/gKVhxxPRR5aqbtBmbNiwQceOHat9+vTRGTNmNH9aKysrv8Dm583CdqolCrDm5+cvqa+vb4gfv3z5ch0/frwuWLBg4zlOP/10nT17dotzN7+Uqv4sDX5OafMAfhL8nP8adizp/PDuhQ5CRMYAL2OLlcep6sxwI4ocRdj0iyadBfX19Tz99NO8+uqrzJs3j2XLllFSUsL3v//9+LSC+7ExTvHNK4fQaORzGNB17dq1dO/e2PY7f/58Bg4cSHZ2Nrm5uVx00UUMHDiQa665htWrV1NdXc3OO+/cWox3YwuoDiAie2B+umuAPgOufbIBM1KKt2XmYWseiW2ZmWKD2WZcdDsQEbkOuBH4EjhAVT8POaQoMRQz1Nns9u2KigrWr18fF8X12HyyR1o7NrBiHFZaWnr5z3/+83Pz8/OzwSYfZGU1Lm9MmzaNL774gsmTJzNu3DhKSko47rjjmp9OgceAs7buv5d5FJWUZX352I3L8gccsGu3/Y/5IKtL4W5seQPSIsy7oQx4oTOIsItuByIi2Vhb0pHBv8dr5hhcdzSClQgOo6Wnw6aoxjKrtgyDPE3NqKZAEubMAMyZM4d7772XWCxGjx49uPPOO1t7fSW2Yt/px9EEW+0nAVdpXc0OZOXkSla71ugVW4Rbh7Vb3r81U7SjgotuByMi/YB3sVHUJap6S7gRRYq9MFHLwcYNbemd/Am2y7CtHAg8p6o9giwYVWXRokXsu+++nHLKKcyYMaPFi2KxWP0zzzwz+7TTTvsF8KZ20o0wganU/wAXkXxTqWnANR013DVMXHRTgIicgPk31ANjVfW1kEOKErlYw/3xwLGYkXkM80lIzIArsYWcP7bn5N98803ftWvXftC3b9/uhYWNmjFlyhR+8pOf0K1b09FylZWVPPfcc5xyyinxb1VjHwzx7cuvq2okzezbQ1FJ2RisjNNh9qmYafrZ5aXFr3TA+UPDRTdFiMivMaezT4EDVTVjb586mJ5YyeEoTIj7YwMQJ2EOYO1CRH7VpUuX/3zooYdqTznllLqsrKwCsIy3WdUBoK62tvbjAQMG3PjZZ5/FHdX2b3ZMPWaAlOgrnDG/62BQwFRSNww2BjwAXJEpgwJcdFOEiORhCwYjsAGGp6v/8ENFRL6Lmd3UZ2dnn1BXV3cINgVhU7fJq4H9SJhyG2zAiJu7j8V2JSZm4IqZqyf6Cn+R5P9KSghrJBYmvJEeiZWIi24KCVpq3sFsCi9V1VZXaJyOR0QOxRbqugCXq+rvg6dOBh7GhDcx1a0ExmFjgTZ33u2AUTSK8KG09AL+gKYinPaz1Xz4a/Jw0U0xInIWNkqoGhipqu+GG1HnQ0T6Y+1ofbFe2580u+vYF5gB9Ev43neAF7fiWgXY3U1chA+jZSa9jKa+wh+m011QUFJ4CbPmDENw41RhI6GOinKpwUU3BETkD1g/6QfA8M6w8JIuiEhXrMwzDBOS47T10fI5WL32a2zab1Le5EGXxEE09RXu0eywz2kqwgvCbDUsKim7k9TVcLdEDGspuzTsQLYWF90QEJFCLNMaCjygqpNCDqlTEHj4PooZnC/FfI+TNTxza2PKxjLrsQmPHZsd9i32QREX4XdUNSXj54MuhWdJD8GNEwOOjWpXg4tuSIjIUKw+WIDNVvtLK4ftA0zA/GHnYTPGnK1ERG4AfoE14Y9U1UXhRtSSYKPGYJqau+/a7LAKzEwpPmFjrqom/XY76MP9kKZllnRhFTAwin28OWEH0FlR1QUi8lNssvDdIvKmqi4Jns7Fms5/hP2OqoJ/J2Jzu5x2EtTSf4E18Z+TjoILGw3xFwePewIRHkDTTHgg1rN8bPCyahF5g6a9whuSEM4UWo5UShd6YdO+fxp2IO3FM90QCd5QDwNnY10No1S1H2aRtwctF1wqsZXxf6cyzqgjIsOx2/N84EpVnRpySNuEiOyM1YLjRj77Njukjpa9wmvac41ga+8qwl042xJVQL+obRl2E/MQCbKaS4CPgWFXX331DExQ96H1XtECbPR386m6ziYItmH/AxOP+4DfhhpQElDVz1T1r6p6qaruh82DOwW4DXgLe18fihnqPwGsFpF3ReR2ETldRJrXjFtjEnZXkM40YAt8kcIz3TRg8ODBo0tKSmadeeaZWV27blFPY9jmiu92fGTRJmjXmoXZO84GjtZOMKlWRLrTtFd4BC17hRfTtFd4efyJYFL2CmxuXbqzCugfJXcyF93w2Qd4sra2tn9ubm7uFo82KrEBjg91XFjRJijdPAScgy1EjlDVr0INKiSCD59Dador3LwboZxAgHc865f1+bsfdIeIdCP9qQBOLS8tfj7sQNqKL6SFhwAXAr8D8nNzc9tT6ikE7sG6H7yjoXX+ExPcCuDkziq4AKoawwz1X4aNW9ITe4VHY6bxRcD3ar9ZQX7RMGjpPZGOFAInAi66zmbJwnZCnccW7PDq6uq466676N+/P3l5eZx44onxpwqwTG4Etr/fCRCR04CbsJ/Luar6fsghpRVBieWN4DEl6BXej0CE83fbb4JkZUVFG7KwuCODlxfC4Szg/9jCgtjjjz/OvffeS35+PuPGjeMvf/kLl112GRMnTowf4t0MzRCRA7HJtIXA1ar663AjihZFJWVZqrohGHWUFH5z5gEcvucO9Oqay4bqeuavXMOUZz5gwWfrknWJGNC1vLQ4EmLm3QvhcCBb2OFzyy23cNZZZ1FVVcX111/P5Zdfzt13382NN97I119/HT9MMOcrBxCRvli7XSHwJ+DWcCOKJHuKSH0yT7hLzwLmfPIN099awbeVNRwxaEf+MPHgZF4iPostEkTlFiLT2IPNfODdcccdPPbYYzz99NP07duX73//+zz88MMMGzaM++67r8k8L9o2mibjCTKzv2P+uq8DP0wn05gIMQjr800a59zbaHM8tF93yi4fw849CsjJEuoakvIrqsPi/igZJ+toXHTD4U7MQrDVeu4333zDQw89xKBBgwAYOHAgH374IQMHDmTkyJHk5eWBlRbOxeu58U6Fe4CRmEn8qR2xLbaTUEBTS8uk8L1RAxi443YctmdvAO6d/XGyBBcs3nTyhtgsXl4Ih39h1oGx1p5ctmwZd911F2vWrOH1118HYOjQoQDk5eVVYy5UR2K30g5cjW2RrgS+E1WT8DQhjw4Q3RP33ZmJIwewZ59urFoTY96ypG4iE8wXORL4Qlp4FAILsZlfTf7I6+vrKS4upnfv3sRiMQ466CAmT54MUCkiL2ACszblEachIvId4HHsZ3iaqv493IiiTVFJ2RnYzr3uyT53l5wsxg7sw93nH0yDKkfe+jIr1rSad7SXdcCF5aXFkfAl8fJCeFQCJ2FDDZuUGbKzs3nwwQdRVT7//HOGDBlCZWUlv/71rxcBE2644Qb/pAREZD/gQUxwJ7vgJoUYSSxZdcnJora+gQaF6roGZi35ig01dXTPz6X/9oXJEl1lE3eN6YiLbri8D1wJ/IZm7WPbb789IsIOO+wQq6mpWT1q1Kie8+fPPxi4+IYbbrgnjGDTCRHpg/kKdMP6lUvDjShjWEISdWFY/57cfs4w3vxkNWtjtRxStD3d83P5uqKa91cm7WYtB4s7EnhNN3zuwSYYNFn4CSbRVgKP5eXlDZo/f/7FwVO3Bxlep0VEumA18QGYGfxF3qmQNJaSRF34Yn01n3y9gdEDd+Cs4f3pUZDLk/9exXnT3mB9ddKaJLKwuCOB13TTg+5Ym9MArNRQg7XB/ACbdACAiNyHbR1eBBySJM/USBF0KtyHuWCtxH4On4UbVWZRVFI2D9smHBXmlZcWDw87iLbimW56sA4zIbkM+DXwe6zZ+9Fmx/0UE9x9gNtTGWAacQUmuDFgggtu8hCRbiJy9oaF/9pOGyJj2tWAdQNFBq/ppg9rgQc2d0CwPfNszOjmByLygqo+nIrg0gERORH7UAL4vqrOCzOeTCAYGX8ScCZwApBf8e+ZFOw5HOmyWVuQdKESeCrsINqDZ7oRQ1XnAz8PvvyDiOwVYjgpQ0SGYFM2soAbVHV6yCFFFhHpLiLfFZHHga+whchTMaP316o+nX+l5HaJSq/zOuDFsINoDy660eQPwGPAdsAjgVVfxiIivbFOhe7AdOC/w40oeohIDxGZKCL/xIT2L9jQ0zzgFeBnQH9VPVzr66ZKVvYULItMZyqBW6NkYA5eXogkqqoicjEwHDgYG9B3ZbhRdQzBB8pjmF/F28AFqhqpN1lYiEhPTFjPxIZYxk3yFauDTgdmqOqqVl5+P/CrFIS5LWSxhZJcOuLdCxFGREZiY2hyMKPuJ0MOKakEnQp3Az8EPsOmP6wIN6r0RkS2p1Foj6ZRaOMLTtOBv7dlAbKopOwOrIMmHX0NYsC08tLiyE0D9kw3wqjqGyIyGRvX/oCIHKCqmeQ6dhkmuNXAKS64rROUX07BhHY8je/rBqzeGRfa9tZprwFOIz1FdzVwbdhBbA2e6UYcEcnCVm+Pw4YwjlfVpPqhhoGIHAs8jd1CfldVfR5cAiKyA7b4dSZwFJAdPFWPbbaZDjyuql9uy3WKSspGAzNJL+GNAceUlxa/GnYgW4OLbgYQjNR+D+iLrez/MuSQtgkRGYx5UvQAfqWq/xVySG2lBzY+qSvwJEn2pQ1+z3GhHUdToX2BRqH9utUTbCVFJWV3YqPO00F4Y8D95aXFl4YdyNbi5YUMQFW/FJHzgeeA60XkZVWdFXZcW4OI9MI6FXpgpuTXhxvRFumF1VAnYX6+MSw7rwAuYht7SEVkJ+wW/0zgCBo7juqAZzCh/YeqfrMt19kCVwAHYLvUkjbGZyuoAt4h4ovGnulmECJyEzAZWAUckOyMp6MRkVxMpI7GMvfRqloRblSt0hMTwklYZluDGe80pwLYlXbacIrIzjQK7VgarT9rsQ/WuNAm1ZR2cxSVlHXDWssGE47wVgEfAKPLS4vT8W+izbjoZhAikoON2T4cKMM6GiLzCxaRO7DFsy8xT4VPQw6pNX6KOZo10LrQJrIB+A/gri2dVET6AadjQjuaRqGtwWqq04F/quqarYo6CQTCO5M2zPhLMjEswz0u6oILLroZh4jsBryL3fZeqapTw42obYjIjzBxqgGOVNXXQg6pNXbDsq32ZHpzsWy4BSKyK41CexiNQlsNPIsJ7ROqmjaG9UUlZV2A27AsPxXCG8N6hq8sLy3OiBFMLroZiIicgtVDa4HDVPWtcCPaPCJyFJZBZWOeCn8KOaRNcSfWt7rJHYCVlZUUFjbxLKjGxPpLABHpD5yBCe2oZsc9jQntk6qatPnkHUHQ1fAo9uHeEeIbA74Fzi4vLX6lA84fGr4NOANR1ceBO7DG+EdEJOmjV5JF4B3xGCa4U9JYcMF+ni0E98033+T8889nr732Yvbs2TQ0deiqW7x48UUi8v9E5A1scOZtmOBWYb7A5wJ9VPVUVX0o3QUXIBDCgcA07P+RrC3DlcH5pgEDM01wwTPdjCUYSf46Vn97BDgv3eq7ItIDeAPYG2uxOiXNe4y/i5VAtgNQVUSEn/3sZ/Tp04fLLruMnj17UldXR05OY2PQG2+8wahRG5PaGFZvnw48laYLhe2iqKSsF9ZSdhXmj1FI+xK6Bkxs1wG3Ag+UlxanbJEw1bjoZjAiMgjzK+iKTVe4L+SQNhIs+j0BHI+NLTpMVdeHG9UW6YGVCTZmu7NmzeKuu+7ikUceAWDt2rX06NGjyYuqq6sZOHDg48uXL38QeDpTzeeLSsqysB1xJ2BdF0MwQa3D6tWC+T4o1q6ahQ1n/RfWtfJi1MxrtgYX3QxHRCYCf8IyrOGqujDkkAAQkduw/s+vMU+FT0IOqa28jPXLAvDxxx9z4okncuutt3LbbbfRs2dPLr74Yo444oiNtd3AB/kK4N5wQg6HQIT3AAZhdd8uWO06hs00W1peWtzpBMhFtxMgIn8EvodllCNUNdTJqSLyA6xmVwscrapRcv7/Xn19/V3Z2dmFAGvWrGHixIkUFBRw4403smTJEh555BGOP/54vve97yW+7nWsQ8Hp5PiOtM7BpdhuqX2BqcCP2vrCIFvZk8ZsJQ9r60rMVtp8SygiY2nsW/1xVARXRPYBzujdu/fZK1euLMzOth243bt3Z4cddmD58uXsvffeDBgwgJUrVzJv3rzmonsQ0BvoyJ1jTgTwTLeTICIHYotWXYCzNjV5IaEuVwyMweaxbakutwizmCwDXtiUCIvI7tj03h2Aqaqa1ts5RWQoje1dQ+Pff+WVV+oOP/zwjQnLO++8w7nnnsvixYsBuOSSSzjppJM4+eSTE09XAVxOBP1fneTiotuJEJFLsaGX64ADE+uowQr0JGwFejts8U1aO88mUGwH1jrgN5gpycYV6KBt7TVMvJ7Bdssl1RBmWwn8e4diInsm9oETZw3wODB9/fr1/bp16zaVhB1pN9xwA8uWLePNN99k77335re//S39+/dvfon7MD8GpxPjotuJCERlBua9+iYwesC1T+ZifrwXYRltMqYRVmIZ8DTgmmW3nFSNCdZJwGJgZLrssgp+JvvRKLSDE55eTSC0wIuqWhN8vxdmqt4l8VwLFiygR48e7Lrrrq1dqhKYAkTaAc7Zdlx0OxnBZIF3gN267nvUwzucdOURdPCuoq+f+M1LGxa89F1MxA5V1Y864FptJhDaA2gU2oEJT3+D7eabDrykqrWbOM0rmMdFW9kq8xsn83DR7YTkdN9hXPdRZ7/Ybb/xkpXbZcsv2EYaaqupmP9Cw7pXHz6hrmL1zA6/YCsEQjsME9kzgMQpyl9jdwDTgZfbWPaYBPyOzZve1GMfPKuBS7CyitPJcdHtZMSdorS+brhk5+Ru8QVJQuvrayU7ey4pdIoKhPZgGoV2j4Snv6RRaP+1FfXlAsxCswdNa9912DbWKmxk/ENYKSfjm/6dtuGi24noDJ6ogdAeQqPQFiU8/QXwN0xoZydhy/FB2MLkgcHXG4AHMaGdiy0uOk4TXHQ7CYEl30vYLXbY7v9vA0cly6ovmBM3gkah3S3h6c9oFNpXO8jbYSjWtbECF1pnC/jmiM7DVCwjC1NwCa4/DHPa2uo5V4HQjqRRaBNbBlbSKLSvqWpH39ov6ODzOxmEZ7qdgKKSsjGYKXY6DBaMEwOObY91XyC0h2FCezqwS8LTKzCLyOnAGykQWsfZKlx0M5yikrJC4EOgX9ixtMIqzDN1k16sIpKNja85AxPanROe/pRGoX3ThdaJAl5eyHymYH246Ugv4BZs7thGAqEdiwntadho+TjlNArt3HTzCHacLeGZbgYTbO1dRfh13M1RBfRbdstJ6zHLxLjQ7phwzCeYyE4H5rnQOlHGM93MZhJp3h+qDfWy/u2yMmyzQp+Ep5bSKLTvuNA6mYJnuhlK4Ba2gqY10LSkbv03rPzfCwD9kEahfc+F1slEfDBl5jKeYJZXupNd0L1mh1Ou+SEwWFUnq+q7LrhOpuKim7kUY/aMaY/k5OZ03Xv0Pi60TmfAa7qZyxja54fbKqWn7sfwAb3o17OAmvoG3l2+hpufWsSHXyZ1F28W1q3gOBmPZ7oZSFDPHZKMc507Yjcqquv453urqKiq48jBO/KnC0fQJSfpfzpDikrKtvlDwnHSHc90M5M9MVvBbea0u17j7U9tAMSuPQt45Zqj2LlHAXvt2I0Fq9Yl4xJxGrC4Q/XadZyOxjPdzGQQZjG4zcQFFyA3yG7rG5Qv1yfFqyaROixux8loXHQzkwKSUM9NpDAvm1vPOACAe2d/zFfJF10hvbwhHKdD8PJCZpJHEkW3V2Eu918wggP79+ShNz/llmcWJ+vUiQjNZo45TibiopuZ1JAkX9ddehbwpwtHsGefbtz58kdMefaDZJy2NRRIevrsOOmGi25mEiNJovu3Hx1G3x75rPi2kvzcbK4/yZoi/vHuSt5bkdQZi4rF7TgZjYtuZrKEJP1u+/Ywr5xdexVy4eG7b/z+wlXrki26OVjcjpPRuOhmJktJ0iJpUUlZMk7TFrKwuB0no/HuhQykvLS4AVgUdhztZGF5abFvA3YyHhfdzGU20RmS2AD8K+wgHCcVuOhmLmXYSPAoUAk8FXYQjpMKXHQzlxeA9WEH0UbWAS+GHYTjpAIX3QwlqOveimWR6UwlcGsQr+NkPC66mc39pP/vOAt4IOwgHCdVpPsb0tkGykuLvwWmkb6bDmLAvUGcjtMpcNHNfK4B0lXUVgPXhh2E46QSF90Mp7y0uBI4m/TLdmPA2UF8jtNpcNHtBJSXFr+C1U3TRXhjwP3lpcWvhh2I46QaF93OwxXAO0BVyHFUBXFcGXIcjhMK4gNYOw9FJWXdgFeAwUB+CCFUAR8Ao8tLi5M62dJxooJnup2IQOhGY5lmqksNMeBtXHCdTo6LbicjELwjsR7eVAlvLLjeUS64TmfHywudmKKSstHAo0AvOmY+WQxrVzs7WMxznE6PZ7qdmEAIB2IbKKpI3pbhyuB804CBLriO04hnug4ARSVlvYALgKuA7kAh7ftQbsDEdh3m+fCA7zRznJa46DpNKCopywLGAycAY4EhmKDWYRN7BfPpVWzySBawEPPDfQp40c1rHGfTuOg6myUQ4T2AQVjdtws2tTeGzTRb6hMfHKftuOg6juOkEF9IcxzHSSEuuo7jOCnERddxHCeFuOg6juOkEBddx3GcFOKi6ziOk0JcdB3HcVKIi67jOE4KcdF1HMdJIS66juM4KcRF13EcJ4W46DqO46QQF13HcZwU4qLrOI6TQlx0HcdxUoiLruM4Tgpx0XUcx0khLrqO4zgpxEXXcRwnhbjoOo7jpJD/D+sFWXzu/wdlAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# n 代表图形 G 的顶点数量\n",
    "n = 4\n",
    "E = [(0, 1, 3), (0, 2, 2), (0, 3, 10), (1, 2, 6), (1, 3, 2), (2, 3, 6)]  # 线段参数(顶点1， 顶点2， 权重(距离))\n",
    "G = nx.Graph()\n",
    "G.add_weighted_edges_from(E)\n",
    "\n",
    "# 将生成的图 G 打印出来\n",
    "pos = nx.spring_layout(G)\n",
    "options = {\n",
    "    \"with_labels\": True,\n",
    "    \"font_weight\": \"bold\",\n",
    "    \"font_color\": \"white\",\n",
    "    \"node_size\": 2000,\n",
    "    \"width\": 2\n",
    "}\n",
    "nx.draw_networkx(G, pos, **options)\n",
    "nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=nx.get_edge_attributes(G,'weight'))\n",
    "ax = plt.gca()\n",
    "ax.margins(0.20)\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 编码哈密顿量\n",
    "\n",
    "量桨中，哈密顿量可以以 ``list`` 的形式输入。这里我们将式（4）中的二进制变量用式（5）替换，从而构建哈密顿量 $H_C$。具体的形式可以通过内置函数 tsp_hamiltonian(G, A, n)直接得到。\n",
    "\n",
    "**注意：** 对于旅行商问题，由于我们总可以选定某一个城市为第一个抵达的城市，故实际所需量子比特数可以从 $n^2$ 降到了 $(n-1)^2$。在我们接下来的实现当中都会使用改进过的哈密顿量来计算。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.237497Z",
     "start_time": "2021-05-17T08:00:16.219567Z"
    }
   },
   "outputs": [],
   "source": [
    "# 以 list 的形式构建哈密顿量 H_C -- 通过内置函数tsp_hamiltonian(G, A, n)\n",
    "A = 20 # 惩罚参数\n",
    "H_C_list = tsp_hamiltonian(G, A, n)\n",
    "# 生成哈密顿量\n",
    "H_C = Hamiltonian(H_C_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算损失函数\n",
    "\n",
    "在最大割问题（[Max-Cut 教程](./MAXCUT_CN.ipynb)）中，我们用 QAOA 构建了我们的参数化量子电路，但除了 QAOA 电路，其他的电路也可以用来解决组合优化问题。对于旅行商问题，我们将使用 $U_3(\\vec{\\theta})$ 和 $\\text{CNOT}$ 门构造的参数化量子电路。这可以通过调用量桨内部的 [`complex entangled layer`](https://qml.baidu.com/api/paddle_quantum.circuit.uansatz.html) 来实现。\n",
    "\n",
    "<img src=\"./figures/tsp-fig-circuit.png\" width=\"900px\" /> \n",
    "<center> 图 1: 旅行商问题使用的参数化电路 </center>\n",
    "\n",
    "上述电路会给出一个输出态 $|\\psi(\\vec{\\theta})\\rangle$，由此输出态，我们可以计算最大割问题的目标函数，也就是旅行商问题的损失函数：\n",
    "\n",
    "$$\n",
    "L(\\psi(\\vec{\\theta})) = \\langle\\psi(\\vec{\\theta})|H_C|\\psi(\\vec{\\theta})\\rangle.\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "然后我们利用经典的优化算法寻找最优参数 $\\vec{\\theta}^*$。下面的代码给出了通过量桨和飞桨搭建的完整网络："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 此处使用内置量子电路：complex_entangled_layer()\n",
    "def cir_TSP(N: int, DEPTH: int) -> Circuit:\n",
    "    cir = Circuit(N)\n",
    "    cir.complex_entangled_layer(depth=DEPTH)\n",
    "    return cir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义损失函数\n",
    "def loss_func(cir: Circuit, H: Hamiltonian) -> paddle.Tensor:\n",
    "    state = cir()\n",
    "    loss = paddle_quantum.loss.ExpecVal(H)\n",
    "    return loss(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 训练量子神经网络\n",
    "\n",
    "定义好了量子神经网络后，我们使用梯度下降的方法来更新其中的参数，使得式（7）的期望值最小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.274144Z",
     "start_time": "2021-05-17T08:00:16.264684Z"
    }
   },
   "outputs": [],
   "source": [
    "DEPTH = 2      # 量子电路的层数\n",
    "ITR = 120      # 训练迭代的次数\n",
    "LR = 0.5       # 基于梯度下降的优化方法的学习率\n",
    "SEED = 1000    #设置随机数种子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里，我们在飞桨中优化上面定义的网络。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.495970Z",
     "start_time": "2021-05-17T08:00:16.496407Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "循环数: 10 损失: 46.0240 用时: 1.7304351329803467\n",
      "循环数: 20 损失: 22.6653 用时: 3.757274627685547\n",
      "循环数: 30 损失: 16.6194 用时: 5.594370603561401\n",
      "循环数: 40 损失: 14.3720 用时: 7.301747560501099\n",
      "循环数: 50 损失: 13.5547 用时: 9.344748735427856\n",
      "循环数: 60 损失: 13.1736 用时: 11.30988073348999\n",
      "循环数: 70 损失: 13.0661 用时: 13.120512962341309\n",
      "循环数: 80 损失: 13.0219 用时: 14.732549905776978\n",
      "循环数: 90 损失: 13.0035 用时: 16.409778356552124\n",
      "循环数: 100 损失: 13.0033 用时: 18.174597024917603\n",
      "循环数: 110 损失: 13.0008 用时: 19.95359969139099\n",
      "循环数: 120 损失: 13.0004 用时: 21.790847778320312\n",
      "得到最小路程: [13.000355]\n"
     ]
    }
   ],
   "source": [
    "# 固定 paddle 随机种子\n",
    "paddle.seed(SEED)\n",
    "# 记录运行时间\n",
    "time_start = time.time()\n",
    "\n",
    "# 总qubit数取为：(城市数-1)**2\n",
    "num_qubits = (len(G) - 1) ** 2\n",
    "\n",
    "# 创建电路\n",
    "cir = cir_TSP(num_qubits, DEPTH)\n",
    "\n",
    "# 使用 Adam 优化器\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=cir.parameters())\n",
    "\n",
    "# 梯度下降循环\n",
    "for itr in range(1, ITR + 1):\n",
    "    # 计算梯度并优化\n",
    "    loss = loss_func(cir, H_C)\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "    # 输出迭代中performance\n",
    "    if itr % 10 == 0:\n",
    "        print(\"循环数:\", itr, \"损失:\", \"%.4f\"% loss.numpy(), \"用时:\", time.time()-time_start)\n",
    "\n",
    "# 显示QNN得到最小路程\n",
    "print('得到最小路程:', loss.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最理想的情况是我们使用的量子神经网络可以找到最短哈密顿回路，同时最后的损失值应该等于这条回路上的权重之和，即旅行商所需要走的最短长度。但如果最后的情况不是这样，读者可以通过调整参数化量子电路的参数值，即 DEPTH，ITR 和 LR，来获得更好的训练效果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 解码量子答案\n",
    "\n",
    "当求得损失函数的最小值以及相对应的一组参数 $\\vec{\\theta}^*$后，我们的任务还没有完成。为了进一步求得旅行商问题的近似解，需要从电路输出的量子态 $|\\psi(\\vec{\\theta})^*\\rangle$ 中解码出经典优化问题的答案。物理上，解码量子态需要对量子态进行测量，然后统计测量结果的概率分布（我们的测量结果是表示旅行商问题答案的比特串）：\n",
    "\n",
    "$$\n",
    "p(z) = |\\langle z|\\psi(\\vec{\\theta})^*\\rangle|^2.\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "通常情况下，某个比特串出现的概率越大，意味着其对应旅行商问题最优解的可能性越大。\n",
    "\n",
    "量桨提供了查看参数化量子电路输出状态的测量结果概率分布的函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.554317Z",
     "start_time": "2021-05-17T08:02:14.500593Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "利用改进后的哈密顿量找到的解的比特串形式： 010001100\n"
     ]
    }
   ],
   "source": [
    "# 模拟重复测量电路输出态 1024 次\n",
    "state = cir()\n",
    "prob_measure = state.measure(shots=1024)\n",
    "reduced_salesman_walk = max(prob_measure, key=prob_measure.get)\n",
    "print(\"利用改进后的哈密顿量找到的解的比特串形式：\", reduced_salesman_walk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因为我们之前为了减少所需要的量子比特数改进了一下旅行商问题对应的哈密顿量，上面显示的比特串缺少了顶点 $n-1$ 的信息，以及所有顶点在时间 $t =n-1$ 的时候的信息。所以我们需要将这些信息加回找到的比特串中。\n",
    "\n",
    "首先为了加上对于 $i \\in [0，n-2]$， $x_{i,n-1} = 0$ 这一信息，我们需要在每 $(n-1)$ 个比特之后加上一个 $0$。接着在比特串的最后，我们为了加上顶点 $n-1$在每个时间的状态，我们加上包含 $n-1$ 个 '0' 的 '00...01'，用来表示对于$t \\in [0,n-2]$来说，$x_{n-1,t} = 0$，同时 $x_{n-1,n-1} = 0$。\n",
    "\n",
    "以下代码通过测量，找到了出现几率最高的比特串，每一个比特都包含了式（1）定义的 $x_{i,t}$ 的信息。我们将找到的比特串映射回经典解，即转化成了 ``dictionary`` 的形式。其中 ``key`` 代表顶点编号，``value`` 代表顶点在哈密顿回路中的顺序，即访问城市的顺序。在以下代码中，我们还将量子电路找到的最优解和暴力算法找到的相比较，从而说明量子算法的正确性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.571954Z",
     "start_time": "2021-05-17T08:02:14.559634Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "参数化量子电路找到的最优解： {0: 1, 1: 2, 2: 0, 3: 3} ，最短距离为： 13\n",
      "经典暴力算法找到的最优解： {0: 0, 1: 1, 2: 3, 3: 2} ，最短距离为： 13\n"
     ]
    }
   ],
   "source": [
    "# 参数化量子电路找到的最优解\n",
    "str_by_vertex = [reduced_salesman_walk[i:i + n - 1] for i in range(0, len(reduced_salesman_walk) + 1, n - 1)]\n",
    "salesman_walk = '0'.join(str_by_vertex) + '0' * (n - 1) + '1'\n",
    "solution = {i:t for i in range(n) for t in range(n) if salesman_walk[i * n + t] == '1'}\n",
    "distance = sum([G[u][v][\"weight\"] if solution[u] == (solution[v] + 1) % n \n",
    "                or solution[v] == (solution[u] + 1) % n else 0\n",
    "                for (u, v) in G.edges])\n",
    "print(\"参数化量子电路找到的最优解：\", solution, \"，最短距离为：\", distance)\n",
    "\n",
    "# 经典暴力算法找到的最优解\n",
    "salesman_walk_brute_force, distance_brute_force = solve_tsp_brute_force(G)\n",
    "solution_brute_force = {i:salesman_walk_brute_force.index(i) for i in range(n)}\n",
    "print(\"经典暴力算法找到的最优解：\", solution_brute_force, \"，最短距离为：\", distance_brute_force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下的代码将字典形式的经典解用图的形式展示了出来：\n",
    "* 顶点中的第一个数字代表城市编号\n",
    "* 顶点中的第二个数字代表旅行商访问此城市的顺序\n",
    "* 红色的边表示找到的最佳路线"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.864346Z",
     "start_time": "2021-05-17T08:02:14.576418Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAADnCAYAAAD7CwxiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABezElEQVR4nO3deViU5foH8O8MM+yCC6CgKOCSC4ob4F6ppWmLmrYeW7U67ac6lS0eq1P6K1tOHc9xqY7HNNPMU5ma+4KZgKK4iwoqigqIgMAwzPL+/ngaZ4YRGWBm3lm+n+vqqst5mXks48v9vvfz3ApJkiQQERERERFRgyjlXgAREREREZEnYjFFRERERETUCCymiIiIiIiIGoHFFBERERERUSOwmCIiIiIiImoEFlNERERERESNwGKKiIiIiIioEVhMERERERERNQKLKSIiIiIiokZgMUVERERERNQILKaIiIiIiIgagcUUERERERFRI7CYIiIiIiIiagQWU0RERERERI3AYoqIiIiIiKgRWEwRERERERE1AospIiIiIiKiRmAxRURERERE1AgspoiIiIiIiBqBxRQREREREVEjqOReABH5JqNRwumSKuQVV6BaZ4TOYITaT4lAtRLxEaHo0DIYSqVC7mUSERG5FPPRs7CYIiKXMBol/HayGJuPFiLzVAmOF1ZAqVBApVRAggRJAhQKQAEF9EYJRklC56hQJMe1xPCuURjcMYLhQUREXof56NkUkiRJci+CiLxXmUaH5bvzsSAtF5VaPapqDGjINx0FgGB/P4QEqDB1aALu6R+L8CC1s5ZLRETkEsxH78BiioicQlNjwMy1R7Bsdz4UCqBaZ2zyewaplTBKwL39YzHttm4I8vdzwEqJiIhch/noXVhMEZHDZeSV4NmlWSjX6FCtb3pI1BaoUiIsSI05D/RFclxLh78/ERGRMzAfvQ+LKSJyGK3egPd+OYwVWWcdcqetPoFqJSb2bYe3b++OABXvwhERkXtiPnovFlNE5BCVWj3+9FU6jpwvd8rdtroEqpToHhOGbx5LRUgAz9QhIiL3wnz0biymiKjJKrV6TJy3E7lFldC6MChMAlRKJESGYMWTgxgYRETkNpiP3o9De4moSbR6A/70VbpsQSHWYERuUSUmf50Ord4gyxqIiIgsMR99A4spImqS9345jCPny2ULChOt3ojDBeV475cjsq6DiIgIYD76ChZTRNRoGXklYjOtzEFhUq03YkVWPjJPlci9FCIi8mHMR9/BYoqIGkVTY8CzS7NccipRQ1TrjHjm2yxoatjOQERErsd89C0spoioUT5YewTlGp3cy7imco0OM39lOwMREbke89G3sJgiogYr0+iwfHe+27Qv1FatN2JZZj7K3DTMiIjIOzEffQ+LKSJqsOW786FQyL2K61MqgO9358u9DCIi8iHMR9/DYoqIGsRolLAgLdftesFr0+iMmJ+WC6ORo/SIiMj5mI++icUUETXIbyeLUanVy70Mu1Ro9diZe0nuZRARkQ9gPvomFlNE1CCbjxaiykNOAtLoDNh8tFDuZRARkQ9gPvomFlNE1CCZp0rgKY0BkgRknuKdNyIicj7mo29Syb0AIvIcRqOE44UVdl0boFJi2m3dcHuvaIQGqHDwXBn+vuYI9uWX2v15dybF4PEh8ejWJgz+KiVW7MnHKyv2N2jNORcrIEkSFO6+I5iIiDxWQ/IRcExGKhTAC8M7497kWLQM8cfJwgp8uP4Yth4rsuvrmY+OwSdTRGS30yVVUNr5TXf67d3xyKA4FFdosf7wBfRt3wLfPJaCFsFquz+va5tmMBglnL5U2dglQ6lQ4PSlqkZ/PRERUX0ako+AYzLyqWEd8eLILtAbJPyy/zw6Robiy8n90Tkq1K6vZz46BospIrJbXnEFVMr6w6JViD8m9YuFwSjhwS/T8fx3+/Bj9jk0C1Tj4YFxdn/eh+uOYcK/dyLtRHGj16xSKpBX3PhijIiIqD725iPgmIz0UyowdWgCAODPS/bg5e+zMS8tFyo/JZ4clmDXezAfHYPFFBHZrVpnhGRHR3iX1s3gr1KioFSDS5U1AIADZ8sAAN2jw5y6xtokANV6z9gQTEREnsnefAQck5HR4YFoGeIPg1HCwYLyRr0H89ExWEwRkd10BiMkO7IiItQfAFBZYz4i1nTCUWSzAKesrS6SJKHGTSfRExGRd7A3HwHHZGRkqLhOozMXQ1V/vJ+978F8dAwWU0RkN7Wf0q7J7sUV4k5biL/5jJuQAD8AQNEVrVPWVheFQgF/Fb/VERGR89ibj4BjMrKoQlwXpPa7+rkhAaoGvQfz0TH4b5CI7BaoVkKB+tPieOEV1OiNiGkedPUOXK92zQEARy6UO3OJNhQAAlV+Lv1MIiLyLfbmI+CYjDxfVo3LVTXwUyrQs214rfe4Ytd7MB8dg8UUEdVPkoD8fMRnp0Ovrf+OV3FFDVZknYWfUoEljw/AF/f1wZ29YlCh1eO/v58GAEzs2w6nZo7FmueG1Pk+t3ZvjdkTe2FopwgAQP+4lpg9sRfu7R9r99L1RgnxESF2X09ERNQg1dWIv3AK+poauy53REYajBIWpOUCAP71QF98PCkJU4fEQ28wYt72k3atg/noGJwzRUTW9Hrg6FFg3z7x19694u8lJegABYwvrwDU9d/JemfVIegNRoztGY24Vq2xN78U7685jJI/Ntua2hL0xrqbzLtHh2FiP3PhFNcqBHGtxDf+Zbvz7frtGCUJHVoF23UtERHRdZWUmPPRlJFHjqCDwfhHPtq3X8kRGTl320kEqv1wT79Y3NErBieLKvDR+mPIuWjfvCvmo2OwmCLyZRUVwP791kXTgQPAtZ4+tWwJZZ8+6OxXg4OoPyy0eiOm/3wI038+dM3Xu7ZpBkCEQV0+23Qcn206bs/vpE5dWodyICERETWMJAGnT9veWDxzxvZapRLKbl3R2VhpVz4CjslIowR8siEHn2zIsesza2M+OgaLKSJfceGCOQxMfx0/jmsePxQfD/TpA/TuLf7q0wdo2xZQKJC86hAO7Txl5wGwdRvUMQI/Z5/DmoMXmvhOdVMogOS4Vk57fyIi8gI6HXDkiG1GlpbaXhsUBCQlmfOxd2+gZ08gONhh+Qg4PyOZj47DYorI2xiNokiqfTft4kXba9VqoEcPc8HUuzfQqxfQvHmdbz+8axSW785HZU3TZlPc9nlak77eHkFqPwzvGuX0zyEiIg9RXg5kZ1tn5KFDwLX2O0VGmrPR9PfOnQG/a7e6OyofAednJPPRcVhMEXkyjQY4eNC6aNq/H6i8xkTzsDDrQOjdG+jeHfD3b9BHDu4YgZAAlUPCwtmaBagwKIF33oiIfI4kAQUFtjcWT9bRNtepk21GRkfD7vPOwXz0VSymiDzFpUvWgbBvnzgownCNb9rt2tm26cXFNSgU6qJUKjB1aAI+3nAM1Tr3HfYXpFZi6tAEKJXsByci8moGA5CTY5uRRUW21/r7A4mJ1hmZlAQ0a9bkZTAffROLKSJ3I0nAqVPWgbB3L3D2rO21SqW5Tc9UNCUlARERTl3iPf1jMXv9Mad+RlMZJWBSA45QJyIiD1BVZT44yfTX/v2iU6O25s1t2/S6dhUt7k7CfPQ9LKaI5FRTAxw+bHs3rfwaQ/uCg82bXk2hkJgoNsO6WHiQGvf2j8Xy3fmo1rvf3bdAlRL3JMciPMh5gUlERE5WWGjbppeTI/YG19ahg22bXvv2DunIaAjmo+9RSNK1jvIiIocrKxObXi2LpkOHxClCtUVFiUCwDIVOnerc9CoHTY0BN87egsIr9Q/xdbXWYQHY+vLNCPJ3n39fRERUB6MRyM21vbFYUGB7rZ+f2O9rmY+9ewMtWrhyxdfFfPQtLKaIHE2SgHPnbNv08vJsr1UoxMlAlm16vXsDbdq4dMmNlXmqBJO/Tner3vBAtRKLH0tF/7iWci+FiIhq02rFjUTLjMzOBq5csb02NNS6YOrTRxRSgYEuXXJjMB99B9v8iJpCrweOHbNtQ7h0yfbagAAxi8KyaOrZ0yGbXuWSHNcSE/u2w4o9Z92inSFQpcTEvrEMCiIid3D5sm0+HjkisrO2mBjbG4sJCWJvsAdiPvoOPpkisldFBXDggHUoHDgAVFfbXtuypW3vdteugMr77l9o9Qbcv2AXDp2+BK1CvraBAJUSPWLCsHTqAASo2L5AROQykgScOWPbpnf6tO21CgVwww22GRnlfTOPrubjuTJoDfL9uM18dC4WU+RVNBoNsrKykJ6ejqCgINx///1ofp0BtHW6eNE2FHJyRGDUFh9v24bQrp3LN73KRpJQ+fKrmFjaHrmt2kGratjcKkcIUCmREBmCFU8OQkiA9xWsRERN5bB81OnEWA7LfNy3TzyFqi0oSAyCt8zHxEQgJKQpvxWPUnkiFxM/Wo/cZlHQqgNc/vnMR+djMUVe5f3338fq1asxcOBAXLp0Ce3atcOMGTOgquuJkNEInDhh24Zw4YLttSqV+Rhy0920pCRx9Kove+cdYMYMVAY3w+S3luJwldKlLQ2BKiW6x4Thm8dSGRRERHVocD4CYh9TdrZ1Rh48KE6irS0iwvYY8s6dvbIjw24XLgBDh6Ly9FlMnvoZDrdsz3z0Qiym6mA0SjhdUoW84gpU64zQGYxQ+ykRqFYiPiIUHVoGc9iZGzIYDPD748S73bt348UXX8TatWvRzLQvqaoK+PZbc9GUnQ1UVtq+UbNmti0I3buLfU9k9umnwEsviZ72ZcugHTce7/1yBCuy8l2y6TZQLXrA3769G1sXiFyE+eiZ6s1HoxHIygLWrzdn5IkT136zjh1tMzImxnc6MuxRUgLceKMoPvv2hXb9RryXdo756IVYTP3BaJTw28libD5aiMxTJTheWAGlQgGVUgEJEiRJfI9QQAG9UYJRktA5KhTJcS0xvGsUBneMYHi4iaqqKixcuBC//vorBg4ciGnTpplf1GhEe4HlH/u2ba0DoU8fIC7OYze9usxXXwFTpoh//s9/gEceufpS5qkSPPNtFso1OqfchQtUKREWpMacB/oimZtpiZyK+eg9rpuPOh3w/vui28BErRZteZYZmZQEhIW5eume5coVYORIICMD6NYN2L5dPLkD89Eb+XwxVabRYfnufCxIy0WlVo+qGgMa8i9EASDY3w8hASpMHZqAe/pzEJpLmVLcgsFgwJ///GdERERg3759+L//+z/07NlTvFhVBbz5priD1qePCIXISBkW7uGWLQPuv1/8+//HP4Dnn7e5RFNjwMxfj2BZZj6UCkDjgDtxQWoljBJwb3Ispo3uxjkZRE7EfPQ+181HANizB1i82HxjsWtXwN/1+2A9mkYDjBkDbN0q9lSnpYmbtpaXMB+9is8WU5oaA2auPYJlu/OhUMAhj1yv/kHuH4tpt/EPssPV1IgjVU192wEBwBtvAOHhdX7J559/jsLCQkyfPh3+DATHWL0aGDdOHG373nvAW29d9/IyjQ7f787H/LRcVGj10OgM1zzHoy6S0QilpEdk82Z4YmgCJvEHMiKnYj56qKIi671Nb74p9vnWgfnoBDU1wIQJIidjYkQhlZBQ5+VNzUcYjTDqtWgRGoTnRnZjPsrEJ3ejZeSV4Nml4hGr1oGPWE13Fpbvzsevhy7wEWtTlJeL/UyWpwUdOmS96bVNG2DGjOu+jVKpRFFRkTNX6lu2bgUmThSF1F//KsK6HuFBakwZmoDHBsfjt5PF2HKsCBl5l2q1CgGSJEGhUEABXG0Vim8ZiMzV30J7ai/2ZG5ASzeacE/kjZiPHsBoFEPga584e+6c9XUDB4oWszpa1pmPDmYwAA89JAqpVq2ADRuuW0gBTcvHLq1DoTt3FNt+mIPbxw7GlKFjXfLbJFs+9WRKqzfgvV8OY0XWWRdu/muHt2/vzs1/dZEkoKDAOhD27gVyc699fefO1nubRo2yCgq9Xo+8vDzs2bMHeXl52LZtG6ZMmYKJEye64Dfj5TIygBEjxLytJ58E/v3vJm02NholnCmpQl5xJar1BtTojfBXKRGo8kN8RAg6tAqGQqHA8OHDsWXLFixevBgPPvigA39DRGTCfHRTWi1w+LDtMeRXrtheGxoqWtdNGTlihNj/+8f3aeajE0kS8MQTwJdfiv1kmzcD/fo1+u3szcc9e/agf//+iI6OxtmzZ6HkXm9Z+EwxVanV409fpePI+XIeSykXgwE4dsz2GPLiYttr/f2Bnj2tTwvq1UucslePH374AQsXLkS/fv0wdOhQDB06lC0MTXXggDiV6PJl4IEHgEWLAD/X/AD0z3/+E8899xwmTJiAH374wSWfSeRLmI9u4vJl22PIDx8WnQC1RUdb31js3VucsFfPD9PMRyeQJODll8XptkFBwLp1wNChLvpoCXFxcThz5gx27tyJgQMHuuRzyZpPFFOVWj0mztuJ3KJKh7Yt2MsnB6ZVVoofwC2LpgMHxMbM2lq0sD1itWtXcYoQye/ECWDIEDHI+M47gRUrXPrf5uzZs4iNjUVQUBCKi4sRHBzsss8m8nbMRxlIEpCfb9umd+qU7bUKBdCli21Gtm7twgXTdf0xaxFqNbBqleiYcaG//OUv+Oyzz/DKK6/go48+culnk+D1xZRWb8B983fh8PlyWYLCJEClRI+YMCydOsD7WhoKC23b9HJycM1dlB062B5DHhvL2RTuKj9f3GE7fRoYPlz0ggcGunwZAwYMQHp6OlauXInx48e7/POJvBHz0QV0OtGRUbtNr6TE9trAQNGBYZmPPXuKcR7knixnLS5fDtx9t8uXkJaWhmHDhiEhIQEnTpyAgj9PuZzX3wZ675fDOCJzUACAVm/E4YJyvPfLEfx9XKKsa2k0oxE4edK2Te/8edtrVSox5NbyblpSkngKRZ6hsFDMyTh9GhgwAPjpJ1kKKQCYMGECiykiB2M+OlhFhW2b3sGDYt9Tba1ambPR9PcuXUR2kmf46itRSAHA11/LUkgBwKBBgxAVFYXc3Fzs378fSUlJsqzDl3n1k6mMvBI89J90l2ymtVegWolvHkt1/1OMqqvF6XmWRVN2tgiL2po1E4WS5ROnHj3E0eXkmS5fBm6+Wfw379VLnOInYyF8/PhxdOnSBeHh4SgsLGSPP1ETMR+b6Px52xuLJ05cuyMjIcG2Ta9tW3ZkeDLLWYuffw4895ysy3nyyScxf/58TJ8+He9YDl0ml/DaYkpTY8CNs7eg8Mo17gjJLKpZALa9crP7zNkoKbE9hvzIkWtvejUNu7VsQ4iPr3fTK3mQigrg1luB338Xd0q3b3eL/vxevXrhwIED+PXXXzHKxT3pRN6E+dgABoMokmq36V28aHutWi1uJFpmZFLSdWchkgeynLX497/bNSLE2datW4fRo0cjMTERBw4ckHs5Psdrnyd/sPYIyjU6uZdxTeUaHWb+egTv3unidgZJAs6csd3fdOaM7bVKpZhPYVk0JSUBUVGuXTO5VnW1CInffwfatxdzMtygkAJEq9+BAwewcuVKFlNETcB8rINGYz44yZSP+/cDVVW214aH256m162bOImWvJflrMVXXwXeeEPuFQEAbr75ZoSHh+PgwYPIyclBly5d5F6ST/HKJ1NlGh1SPtgoex/49QSolMh4Y6TzJlXrdOLpUu3TgkpLba8NChKtXJZ303r2BHhqmm/R6YBJk8TeqNatxeT2zp3lXtVVpl7wqKgoFBQUwM9FR7MTeRPm4x+Ki23b9I4eFXuDa4uNtW3Ts5jfRD7CctbiU08B//qXW/0ZmDx5MhYvXoxZs2bhtddek3s5PsUrn0wt353vTn++r0mpAL7fnY8pQ68/Hdsu5eXi7pll0XTwIFBTY3ttZKRtm17nzi6bGURuymgEHn1UFFItWgDr17tVIQUAPXv2RMeOHXHy5Ens3LkTQ100x4PIm/hcPkoSkJdne2Px7Fnba/38bNv0evcWh0WQbztwABg9WhRSDz4IzJnjVoUUILo3Fi9ejJUrV7KYcjGvezJlNEoYMGuTW/aC1xbVLAC7Xh8BpdLO/yElSWx6rd2md/Lkta/v1Mm2DSE62u2+AZDMJAl4+mlg7lxxBO+mTUBqqtyruqZXX30VH330EV588UV8+umnci+HyKN4dT4C4gbi4cO2+5vKy22vDQkRreuWGdmjh+jUILJkOWvxrruA7793yzmYVVVViIiIgEajwZkzZxAbGyv3knyG1z2Z+u1kMSq11zg4wQ1VaPXYmXsJQzpF2L5oMIhZTbXbEIqKbK/19wcSE62Lpl69gLAwp66fvIAkAa+/LgqpgAAxcNBNCylA3Hn76KOPsHLlSnzyySecp0HUAF6Tj4BoWa99DPnhw6JdubY2bWxvLHbsyI4Mql9+vhgRcvGiaPH77ju3LKQAIDg4GLfddhtWrlyJH3/8Ec/JfMKgL/G6Ymrz0UJU1RjkXoZdNDoDNh8txJCYYPEI2fJu2v79YjNsbc2b2/Zud+vmtv9zk5ubORP48EMx22TFCnEcuhtLSUlBTEwMzpw5g6ysLPTr10/uJRF5DI/Mx46tREte7Ta9vDzbL1IoxAmktdv02rRx5dLJW1jOWhw4EPjxR9lmLdprwoQJWLlyJVauXMliyoW8rpjKPFUCT+lblCQg85ftwF09r73ptX172/1N7duzTY8c45//FEe6KhTAN98At98u94rqpVQqMX78eMyZMwcrV65kMUXUAB6Xj+t3AY8MAy5dsr0gMFAclGSZjz17AqGhrl4qeaPLl8WIkJwc0Q66erVH/NkaO3Ys1Go1tm/fjqKiIkRGRsq9JJ/gVcWU0SjheOE1hspew2OD4zCpXyy6tG4GP6UCn23MwWebjjfo816+pQtGdI1Cu5bi1Lsj58vx0bpj2H36st3vkRMaBUmhgMIUCqbiKSkJaOkBgwvJM/33v+Yhg/PmAffdJ+96GmDChAmYN28eCgoK5F4KkcdoSD4CjsnIO5Ni8PiQeHRrEwZ/lRIr9uTjlRX77f76HP/mkC5dgqJlS3M2mv5+ww3iiTqRo1VUAGPHijbSLl3EgUwyDq1viObNm2PEiBHYu3cvcnJyWEy5iFd9JzpdUgWlnU9tEtuGo0yjw/kyDdq1aNwR4OP7tMWVaj3WHriApNhwpMa3wn8eScaIT7bZvcFXGRiA02cKERfDwolc5IcfgMceE//88cfA1KnyrqeBbrrpJpSWlkKv10OSJO6bIrJDQ/IRcExGdm3TDAajhNOXKtG5dbMGf70ywB+nDxxHXI+O7Mgg16g9a3HjRo+br7lkyRIEBwfDy86Xc2teVUzlFVdAZefJPy8tzwYAzP9Tv0YHxfPf7UPWGfEUKtjfD5lvjESzQDX6tG+BdYcu2PUeKpUKeeV6xMU0aglEDbNuHXD//aKtdPp04KWX5F5RgymVSoSEhMi9DCKP0pB8BByTkR+uOwYAmH5790YVUyq1GnmqMMSxkCJX0OlEl8amTWLW4saNYsaYh2nJriaXU8q9AEeq1hkhubAj3FRIAYACgMpPfMO/UHaNgyPqIAGo1nvGhmDycDt2AOPHi8B48UVgxgy5V0RELuLqfHQE5iO5TO1Zixs2uN2sRXJfXlVM6QxGyPFU00+pwOxJSQhQ+eGX/QXIPltm99dKkoQaN55ET14iK0v0gGs0osXvk0+8tm2mpKQEu3btknsZRG5FrnxsCuYjuYQkAc88AyxZIg6ZWLtWHGbipZiRjudVxZTaT+nynw8D1UosmNwftyVGY9PRi1dbI+ylUCjgr/Kq/wzkbg4fFqcSlZcDkyYB8+d7bSEFANnZ2XjkkUfkXgaRW5EjH5uK+UhO52GzFh2BGel4XvVdKlCthAKuS4vwIDWWPD4Aw7tG4Yess3jimz2oMTTsLpoCQKCKgwPJSfLygFtuEUcL33YbsHix1wyqPHnypM2v5ebm4sKFCzhx4gR+++03GVZF5J5cnY+OwHwkp6s9a/Gmm+RekcMwI13Hqw6giI8Ihd5oXx/Dvf1jkRzXAj3ahgMAbu3eGu1aBGH94YtYf/giJvZth9mTknC4oAxjvthxzff46qH+6NehBUqralCu0eHNMd0AANtyirAtp8iudeiNEuIjuJmenODcOTGxvaAAuPFGcYqfv7/cq3KYp556Cq+99hp69uyJ7777Dps2bUJeXh7i4uIwZ84c9OnTR+4lErmNhuQj4JiMvLV7a9zavTWS2jUHAPSPa4nZE3sh89RlLNudX+8amI/kVJazFhcv9ohZiw3BjHQdryqmOrQMhtHOpvDkuBaY2M98Skv3mHB0jwnH2csarD988Wo7xPXCp024mITdPNgfjw6Ov/rr5Rqd3cWUsboaHf5vhnisnJIiTo7xtF4Mcj/FxeKJVF4ekJwM/PwzEBQk96oc6uabb8Ydd9yBXr16ISYmBuPHj8e4cePQwkPmgRC5UkPyEXBMRnaPDrN6j7hWIYhrJYoje4opo0aDDs88JvIxNVXMmApu3MmCRFYsZy3Onw/ce6+863ECZqTrKCQvO4j+9i/ScLCgvMnv8/bYbnh8SAKeXrIHaw7ad8x5Y/Q8fxyr/vsX8y+0aWMurFJTgf79gfBwp30+eaGyMvFEas8eoEcPYNs2oFUruVflcHl5eejevTs0GuvTM48ePYqjR49i2LBhPCKWyIKj8hFwTUba5KOfH9Crl3VGdu0KKL1qxwI528qVYv+w0ShmLXrgiBB7MCNdx6ueTAFAclxLHCoob/IBsIM6RuDn7HNOLaQUAJL7dQbavw1kZIi/LlwQR3P+9NMfFylEWFiGR8+egFrttHWRB6uqAu64QxRSHTuK4129sJACgPj4eDz++OO4dOkSiouL8fnnn2PVqlXw8/PD4MGDsWLFCvTu3RuvvPKK3EslcguOykfA+RmpUADJN/UBhiwQ2ZieDhw8COzdK/6aO1dc2KyZePpumZHR0U5ZE3mBdevELCmjEfjb37y2kAKYka7kdU+m0o4X4anFe1BZ4/6zKYL9/TB/cn8M6RQhfkGSgBMnRGikp4sA2btXzAWyFBgI9O1rbn1ISQHi4tge6OtqaoC77gJ+/RVo21bMlYqLk3tVTqXVarF371688847iIuLw8MPP4wBAwYAANLT0zF16lTs379f5lUSuQePzkcAqKgQYx4sMzL/Gu2CsbHmwio1VeRlaKjrFk/uaccOcbKtRiNmLXrxiBATZqRreN2TqcEdIxASoPKIsGgWoMKgBIunBgqFGBLXuTPwpz+JX9NqgX37zHfmMjKA48eBnTvFXyaRkdZ35pKTxeA58g16PfDAA6KQiowUk9u9vJACAJ1Oh3/84x946KGHcP/991u9FhwcjJCQEBQUFCAmJkamFRK5D4/OR0AURMOGib9Mzp83Z2N6OpCZKQqs/Hxx6A4g2gATE60LrO7dveZkU7KD5azFxx/3iUIKYEa6itc9mQKABWm5+HjDMVTr3HfYX5BaiZdvuQFThiY0/ItLSsxtgaY7dJcu2V7XpYt1gZWU5FWnudEfjEYRDgsXiv11W7aIjdo+IikpCUuWLEFiYuLVX1u1ahW+++47TJw4EePHj4ckSVD4QHAS1cfr89FoBI4etS6w9u8HDLUKyJAQsSfZMiPbtXPMb4Dcy5EjogAvLgbuuQf49lufKqSZkc7nlcVUmUaHlA82QuvGk9MDVEpkvDES4UEO2PskSeLUNsvWh6ws8VTLkr+/+CHbsj2wY0efuDvjtSRJtCt8/rk45WrDBmDQILlX5VLvvPMOMjIykJqaiqysLOzYsQPt2rXDU089hYceegjBPP2L6Cqfy0dA7CXdu9c6I0+dsr0uJsZcWKWkiA6PZs0cswaSR14eMGSIGBEyZgzwv//53E1lZqTzeWUxBQDTfzqI5bvzUe2GgRGoUuKe5Fi8e2di/Rc3Vk0NcOCAdXgcPWp7XcuW1q0PKSlee2CBV3r7beDvfxfh8Msv4jh0H1NdXY38/Hz8+OOPCA0Nxfjx49GmTRu5l0Xktnw+HwHg4kXREmiZkWVl1tcoFKId0DIjExPFgFdyfwUFwNChQG6umLW4dq3XjQixBzPS+by2mNLUGHDj7C0ovKKt/2IXax0WgK0v34wgfxc/Zi4tFeFh2R5YWGh7XceO1q0PvXuLQy/IvXz0EfDqq6JdYcUKYNw4uVfkNoxG8UOikkcmE9lgPl6D0Sj2I1u2B2Zn2x4AFRQE9OtnnZHt27PDw90UF4sC6vBh8YRx40YgLEzuVbkNZqRjeW0xBQCZp0ow+et0t+oND1QrsfixVPSPc4Oz/SUJOHPG+s7cnj1ig6YltVrst7J8etW5M2d7yGnePOCpp8Q/f/ON+cASH8e+byL7MB/tUF0tDoCyzMiTJ22va93atj2weXNXr5ZMLGctJiYCW7ey4+YPzEjn8OpiCgDe+vEAVuw56xbtDIEqJSb2i8Xfxzm5faEpdDoxy8Py9MDDh0XhZal5c/NsD1OAREXJsmSf8+23oniSJGDOHODpp+VekVsxGAzw86HNxUSNxXxshOJi2/bAkhLb67p2tW4P7NnT5/bqyKKqChg9GkhLE102aWmcO1aL0WjkEykH8/piSqs34P4Fu3CooFzWDbcBKiV6xIRh6dQBCFB52A965eXA7t3W7YHnz9teFxdn3frQt69P9ic71c8/AxMmiJOpZs0CXntN7hW5DYPBgJKSEvz000945JFHoOK+BqLrYj46gCSJp1WW7YF794p9y5YCAkQmWhZY8fFsD3Qky1mL7dqJQsoHRoTYS6fTISsrC5cvX8bo0aPlXo5X8fpiCgAqtXpMnLcTuUWVsgRGgEqJhMgQrHhyEEICvOAHPEkCzp2zvjO3ezdQWWl9nZ8f0KuXdYHVtSvbAxtr0yYxJ0OrBaZNAz74QO4VuRVJktC1a1fk5ORgy5YtuOmmm+ReEpHbkzsf1UqgU+tm3pOPgPgevX+/dYGVk2N7XUSEdXtgSoo4FIoaTq8H7rtPzBaLjAS2bxc/b9BVa9euxZgxY9CvXz/s3r1b7uV4FZ8opgARGJO/TsfhgnKXtjQYdVpEB+ixefoE7wmKa9HrxSwHywLr4EGxqddSs2a27YF8BF+/338XJ/VVVgLPPAN88QXvaF7DtGnTMGvWLDz33HP4/PPP5V4OkUeQMx9RcgYb3hqHGxLiXPa5srh82bo9MD1dtAzW1rmzdYHVu7d4qkV1qz1rcetW8e+NrGg0GkRGRqKyshKnTp1Chw4d5F6S1/CZYgoQLQ3v/XIEK7LyXbLpVq2UULJ7DUo2LcAPy5dhwoQJTv9Mt1JRITaAWrYHnj1re11srPXTq379xEBFErKzgZtuEqcxPvQQ8J//8OleHTIzM5GSkoK2bdvizJkz7AsnspOr8zFQpYT/uSwcWDQDSYk9sGPHDoSGhjr9c92GJIlZV5b5mJUlDr2w5O8vCgPL9sBOnXgzzYSzFhvknnvuwffff49PP/0UL774otzL8Ro+VUyZZJ4qwTPfZqFco3PKXbhAlRJhQWrMeaAvNi9bgNdffx3BwcHYsWMH+vTp4/DP8ygFBSI8TAGSmQlcuWJ9jVIpTuCxLLC6d/epieVX5eSIORmFhcD48cDy5Zxxch2SJKF9+/Y4e/Ys0tPTkZKSIveSiDyKK/OxYxgwYMAAHD9+HOPGjcMPP/zg2zdAdDrzfEhTRh45Yntdixa27YGRka5frzuwnLW4ejUwcqTcK3Jr3333He6//34MHToU27dvl3s5XsMniylAzNmY+esRLMvMh1IBaBxwJy5IrYRRAu5NjsW00d0Q5O8HSZLw6KOP4r///S/atWuHjIwMRLOtzcxgAI4ds24P3L9f/LqlkBCgf3/rAqtdO3nW7CqnT4tCKj8fuPVWcfgE2z3q9cILL+Dzzz/Ha6+9hlmzZsm9HCKP46p8BIBjx45hwIABKC0txbRp0/AB94JaKysTe5ItC6wLF2yvi4+3bp/v08f7D4DirMUGKy8vR2RkJHQ6Hc6fP4/WrVvLvSSv4LPFlEmZRofvd+djflouKrR6aHQGm1PAr0ehAILUfggNUOGJoQmY1D8W4UFqq2u0Wi1GjBiB3377DSkpKdi6dSuCvP2bXFNUVYl2B8vj2U+dsr0uJsa69aF/f7EnyxtcuCAKqRMngMGDgXXr2Ppop23btuGmm25Cp06dkJOTw5kaRI3kinwEgI0bN2L06NEwGAxYtGgRJk+e7MDfhZeRJHGDzbI9cM8ekZuWVCoxH9IyI7t08Z4Wcc5abLTbb78dq1evxrx58/DEE0/IvRyv4PPFlInRKOG3k8XYcqwIGXmXcLywAkqFAiqlAhLMg84UAPRGCUZJQpfWoUiOa4XhXaMwKKEVlMq6f2grKipCSkoKTp06hfvuuw/ffvstf8hriIsXrdsDMzLEHTtLCoVoB7R8epWY2KC2uOPHj2Pjxo3o3bs3+vXrB3855oKUlIg9UgcOiLuLW7aITbVkF4PBgOjoaBQVFeHAgQNITHTzuTVEbs7Z+QgA//73v/H000/D398fW7ZswSDue7GfXg8cOmT99OrQIdv5kOHh4gAoywKrgU8m3CIjLWct/utfwJ//7Po1eLCvv/4ajz/+OEaNGoVff/1V7uV4BRZTdTAaJZwpqUJecSWq9QbU6I3wVykRqPJDfEQIOrQKbnAxdPDgQQwcOBAVFRV499138fbbbztp9T7AaASOH7duD8zOFj3nloKCxIEWlgVW+/bX3LxrNBqxdetWLFy4EPn5+VAoFPj4449du8/tyhXR852RIY513b7dd3vhm2Dq1Kn48ssv8c4772D69OlyL4fIqzgjHwHgueeewz//+U9ERkYiIyMDcZwR1HhXrognVpYF1rlztte1b2/dHtivnzjI4RrcIiM5a7HJiouL0bp1ayiVShQVFaF58+ZyL8njsZhysdWrV+OOO+6AJElYvnw5Jk2aJPeSvEd1NbBvn3WBdfKk7XWtW4vQWLHCZiK9VqtFwB/7kv7617/Cz8/PdftuJAl44QVx7HlcHLBjB9C2rWs+28uY5mkkJSVh3759ci+HiOyg1+sxduxYrF+/HomJidi5cyeaeUvrtjs4d866PXD3bnHqriU/P6BnT5FFDz4IqG23LciWkZcuAQkJQHk58MYbwPvvu+ZzvdDw4cOxZcsWfPPNN/gTWySbzEuaZz3H2LFjMXv2bADAww8/zMFpjhQYCAwYIELg22/FfqOiInHCz9/+BoweLQYiXrwoZmDVfooFICAgANV/HE2r1+vRqlUrGGvPynIWhQKYOVMMHty4kYVUEwwfPhxhYWHIzs7GyWsV1ETkdlQqFZYtW4YbbrgBBw8exAMPPABD7cOIqPHathWnws6aJdrHS0tFO/lXXwFPPCH2WEmSuClpNNoUUoDMGRkUBKxdC7z0kjjBjxrNNKpn5cqVMq/EO/DJlAwkScLUqVPx1VdfITo6GpmZmWjLH5xdQ5LE06oLF8SBFYGBVi8bjUYolUrs3r0bd999N9auXYvu3bubv9YB+9xOnjwJhUKBhISEa19gNHrPJmEZPfjgg/j222/x0Ucf4ZVXXpF7OURkpxMnTiAlJQWXL1/GK6+8go8++kjuJfmOykpxAFRiojiCvZbrZqSDXDcja2pEkcc9501y9uxZxMbGIigoCEVFRQjhAVdNwp/YZKBQKPCvf/0LN954I86fP48777wTVbVP4iHnUCjEwMMhQ2wKKUAcXjBnzhw89dRTmD9/vnVIVFcDjzwCvPuuOF2vpKTBHz9v3jw8+uijuOuuuzBjxgxcqT1jC2Ah5SC880bkmTp16oSVK1dCpVJh9uzZ+Prrr+Veku8ICREnyV6jkALqyci0NGDKFGDBArGHWa9v8MfXm5H+/iykHKBdu3ZITU2FRqPBunXr5F6Ox+OTKRldunQJqampOHnyJCZOnIhly5b59sBCmeXk5GDhwoUoLS3Fiy++iC5dulw9pQqAaImoHTCdO1ufjJSUVOcsqJ9//hlvvvkmNm7ciKqqKjz99NMYO3Ysnn32Wef+xnxUZWUlIiIiUF1djXPnziEmJkbuJRFRA3z55ZeYOnUq1Go1Nm7ciGHDhsm9JJ923YysrgbefBP45BPzFwQHiw4Qy4xs167OYogZ6VoffvghXnvtNTz44INYvHix3MvxaPzJXUatWrXCqlWrEBYWhhUrVmDGjBlyL8lnXblyBffddx8WLFiAoUOHXj1FyupEqoAAYOlS4C9/AQYNEk+2jh8HliwBnn9eBEVYmPj7888DWq3VZ1y4cAHTpk1D69atER8fj5deeglHrjXdnhwiJCQEo0ePBgD8+OOP8i6GiBpsypQp+Mtf/gKdTocJEyYgNzdX7iX5rHozMiBAdG588onY9xsfL2Zfbd8OzJ4NTJokTg6MiRHDdX/80ebJFTPStcaPHw8AWLVqFWpqamRejWfjkyk3sG7dOowZMwZGoxFLlizBAw88IPeSfI7pyNc9e/YgIyMDBw8exEsvvYQpU6bUfcSvTic271qeHmj5jX/VKnHM+R/thJWVlTAYDAgLC4PRaERaWhpmzJiBLVu2AADOnz+P6OhoZ/9Wfco333yDhx56CCNGjMDGjRsBiGOdT5dUIa+4AtU6I3QGI9R+SgSqlYiPCEWHlsH1zsQhItcwGAy48847sWbNGnTr1g2///47wjl3z+UalZGFhUBmpnVGlpaK17p0EXuzLPbqMCNdr1evXjhw4ADWrl2L0aNHMx8bicWUm/jiiy/w/PPPIyAgANu2bUNqaqrcS/J5VVVVCK5j3kadyspEeGRkiELr/ffFnbhr7M/Kzc3FCy+8gFWrVmHGjBnQ6/X4O08ocqjLly8jKqo1/Nv3xPOz5uPAhapaA0elq+eKKKC4OnC0c1QokuNaYnjXKAzuGMHwIJJReXk5Bg4ciMOHD2PUqFH45ZdfoGrAMHZyjgZnpCSJbg7T8exhYaI1sI73YEY639/+NgMfLvoJyXc9irCOfZiPjcRiyk1IkoSnn34ac+fORevWrZGZmYnY2Fi5l0VNJUniCZZSCdQK/+rqakydOhVRUVHIyMjA6tWrERYWJtNCvU+ZRoflu/Mx68dM6CQllAFBAOz/pq8AEOzvh5AAFaYOTcA9/WMRHmR7VDAROV9ubi5SU1NRXFyMF154AZ999pncSyJHMOWjn5/NS8xI5zHl4782H8Olsgoo1IFQNGDPPvPRGospN6LT6TB69Ghs3rwZSUlJ2LFjB0JDQ+VeFjmKXn+1oJIkCWVlZejUqROio6OxZcsWREREWF9fXQ2cPi32ZPXvL/ZitW4tw8I9i6bGgJlrj2DZ7nwoFEC1rukzUILUShgl4N7+sZh2WzcE+dsGPxE5V1paGkaMGAGdToe5c+fiySeflHtJ5Ci1Ro/YlZEaDTB/vsjVlBRxAJS/v4sX7lmYj87BYsrNlJSUYMCAATh+/DjuuusurFy5kif8eYPt28UduAEDrHrE582bh1GjRl3dzHtVTY2YhzVggJj2btKhg/XJSH371tki4Ysy8krw7NIslGt0qNY7fpBkoEqJsCA15jzQF8lxLR3+/kR0fQsXLsSjjz4KlUqFdevWYfjw4XIviZrKYADeew94+WWgWTOrl+rMyMpKcf28eeZfCwgA+vQxZ2RKCtCxI49S/wPz0XlYTLmhnJwcpKamorS0FK+//jpmzpwp95KoKTIzgeHDxZOmjAygd++r39ytjl6v7fx54Pffzf3lu3cDFRXW1/j5AT17moMjNRXo2vWaLRPeTKs34L1fDmNF1lmH3GmrT6BaiYl92+Ht27sjQOVb/66J5Pbqq6/io48+QosWLZCeno7OnTvLvSRqLEkCnngC+PJLMf9xyxarlvg6M1KvB/btM+djejpw7Jjtda1aiWy0LLBatXLe78cNMR+dj8WUm9q0aRNGjRoFg8GA//73v3jooYfkXhI1xsGDwI03igG/DzwALFrU+ELHYBCnBZpORUpPF4dcGGt9c2zWzNwWaAoQL56xVKnV409fpePI+XKn3G2rS6BKie4xYfjmsVSEBHAzPJGrGAwGTJgwAT///DO6dOmCXbt2oUUdQ2bJjUkS8Mor4jj1oCBg3ToxMLixSkvNpweaMrKw0Pa6jh3N3R0pKeIG5zUOifIGzEfXYDHlxubOnYs///nP8Pf3x+bNmzF48GC5l0QNceKECIYLF4A77gB++AFQO3iDZmWlOF7W8ujZM2dsr2vXzvrOXP/+gBfsx6vU6jFx3k7kFlVC68KgMAlQKZEQGYIVTw7yicAgchdXrlzBkCFDsH//fowcORJr1qyB2tHfX8m53n0X+NvfRC7+/DPwx1xAh5Ekse/YVFhlZAB79oi9VpbUarHfyrLA6txZHIzhwZiPrsNiys09//zz+OKLLxAZGYmMjAzbvmFyT2fPipaF06dFi9/q1a6783X+vAgNU4BkZlrvuwJESPToYb3/qnt3mxMH3ZlWb8B983fh8PlyWYLCJEClRI+YMCydOsBnWhqI3MHp06eRkpKCwsJCPP3005gzZ47cSyJ7ffYZ8Je/iCxatgyYONE1n6vTiY4Ry/bAI0dE4WWpeXPb9sCoKNes0QGYj67FYsrN6fV6jB07FuvXr0diYiJ27tyJZrU2aJKbKSwEhg0T/dsDBgAbNsj7FMhoFGuxbH3Yv99m+jxCQoB+/azbA9u1c9vNu2/9eAAr9px1aetCXQJVSkzsF4u/j0uUeylEPuX333/HTTfdhJqaGvzzn//EM888I/eSqD5ffw08/rj45//8B3jkEVmXg/JysSfZMiPPn7e9Li7OOh/79hXtiW6I+ehaLKY8QGlpKQYOHIijR49i7Nix+Omnn+DnYwcMeIzSUuDmm8XG2F69gK1bAXfs5ddogL17rcMjL8/2ujZtrFsfkpPFoEWZZeSV4KH/pLtkM629AtVKfPNYqs+dYkQkt8WLF2Py5Mnw8/PDmjVrcOutt8q9JKrL8uXA/feLm3z/+Afw/PNyr8iWJAHnzlm3z+/eLdrqLalU5gOgTBnZtavs7YHMR9djMeUhTpw4gdTUVJSUlODll1/G7Nmz5V4S1VZZCdx6K7Bzp+i3TkvzrLlQhYXmzbumACkttb5GoQC6dbNuD0xMdPxesOvQ1Bhw4+wtKLyiddln2iuqWQC2vXKzT87ZIJLTm2++iQ8++ADh4eHYtWsXunbtKveSqLY1a4C77hJdEe+9B7z1ltwrsp9eDxw+bN0eeOiQ7QFQYWHmA6BMBVZ0tMuWyXyUB4spD7J161bccsst0Ov1+PLLL/G46TE5ya+6WhwysXEjEBsL7NgBtG8v96qaRpKA48etn17t2yd6zi0FBYl2B8sCq0MHp7UHvv3TQXy/O98t2hdqC1QpcU9yLN6903vbGYjckdFoxMSJE/G///0PHTt2RHp6Olr52BHYbm3bNnHARHW1OMHvww/dtoXcbhUV4kALywLr7Fnb62JjrdsD+/WzmjfpSMxHebCY8jBffvklpk6dCrVajQ0bNuDGG2+Ue0mk1wOTJgE//iieRKWliSdT3kirFQWVZYF14oTtdVFR1ht3U1LEht4mKtPokPLBRlk31NYnQKVExhsjER7Ek8WIXKmyshJDhw7F3r17cdNNN2HdunXw9/eXe1lkmrVYUSFmSs2d6/mFVF0KCqxPD8zMBK5csb5GqRQdHZZPr7p3b/J8SOajfFhMeaCXX34Zn3zyCVq1aoXDhw8jyoNOmPE6RiPw8MPA4sWiWNi2TeyV8iWXLtm2B166ZHvdDTdYP73q1Qto4A86C9Jy8fGGY27VC15bkFqJl2+5AVOGJsi9FCKfc/bsWSQnJ+PChQt44YUX8Omnn9Y9GJ2cz5GzFj2RwQAcPWr99OrAAfHrlkJCbNsD27Vr0EcxH+XDYsoDGQwGjBs3DiNHjsQTTzyBIDc9TcbrSRLwzDPAv/8tvhFu2iS+Cfo6SQJyc62fXu3dK55qWQoIAPr0sS6wEhLqvGNpNEoYMGuTW/aC1xbVLAC7Xh8BpZI/xBG5WkZGBqZMmYK1a9eidevWUHnQyAev4opZi56oqkrMh7QssE6ftr0uJsa6PbB/f6CO05yZj/JiMeWh9Ho9ampqEBwcLPdSfNe0acCsWaIoWLtWnOJH11ZTI45jtyywjh2zva5VK9v2wD/2PaQdL8JTi/egssZg+3VuJtjfD/Mn98eQThFyL4XIJ+l0OiiVSp58Kxc5Zy16oosXrdsDMzKAsjLraxQK0Q5oWWAlJgIqFfNRZiymvFRZWRnCw8Oh1+t5V84ZZs4E3nhDHI26cqW460YNc/myebaH6a+iItvrOnUCUlLwTtcxWFjVHJ7wDUuhAB4dFI/pt3eXeylEdA3MSCdyt1mLnshoNB8AZSqw9u2znQ8ZFAT064d3Bj6IhX7tIMH9n/Z4Yz6ymPIyGo0GS5cuxeLFi5GYmIjy8nIsXLhQ7mV5lzlzgGefFd8Rvv0WuO8+uVfkHSRJ3MW0fHq1Z484/QnA7Y98hoNtOsm8SPv1bBuGVc8OlXsZRGSBGelknjJr0RNVV4uWecv2wNxcAMxHubGY8iIlJSWYMWMGfv/9d3zxxRfo2LEjXn75ZVRXV2P58uVyL887LFokDpwAgPnzgalT5V2Pt9PpgIMHYdyVjm6nY6BV1N+y8/GkJAzuGIEWIWpUag04cK4UH/56DIfOl9v9sa+OugF3JsUgMjQA1Xojjl24gk835uD33GscrFGHAJUSR98dzc3vRG6CGelknj5r0RMVF8OYno5uO4zQwr5hwY7IyJdv6YIRXaPQrqXYanLkfDk+WncMu09ftuvrvS0f5R3TTA712GOP4fvvv0diYiKysrIQGRmJRYsWoWPHjrh48aLcy/N8//sf8Oij4p9nz2Yh5QpqNdCnD05PmgylnSf/tW0ehPS8S/h+91lcrqrBjV2iMG9yvwZ9bGyLYGSfLcPyPWdxpqQKKfEt8fXDyQhS27//QqlQ4PSlqgZ9LhE5DzPSibRaYPx4UUjFxoqZiyyknC8iAqdTb4KyAQd7OCIjx/dpC4VCgbUHLqCgVIPU+Fb4zyPJiGoWYNfXe1s+slHYS7z77rsICgrCuXPnoFQqMWTIEMTHx+O2227D66+/jrCwMLmX6NnWrxftfEYjMH068PLLcq/Ip+QVV0Bl58k/9y3YdfWfe8SEYfVzQxEdHgSVUgG90b4H8c99t/fqP4cHqZE9/VYE+fshItQf+Zc1dr2HSqlAXnEl4iKcM5yRiOzHjHQivV7k44YNooDatMnzh9Z7kIbkI+CYjHz+u33IOiOeQgX7+yHzjZFoFqhGn/YtsO7QhXq/3tvykcWUl6iursZrr70GpVKJ8vJydOzY8er09/DwcJlX5+F++w0YN06cSPfCC8CMGXKvyOdU64yQGnD0xEMDO6BzVDMM6ij+H1iQlmt3SJjcmRSDfh1aoG970e//y/4CuwspAJAAVOvd/2QlIl/AjHQSoxF47DExtL55c3Hj0VuH1ruphuYj0PSMNBVSAKAAoPITxdyFMvsy0tvykcWUl2jVqhVeeuklvPXWWygoKMCVK1fQrI55BNQAWVnAmDGARiNa/D75xHsnt7sxncGIhuzuHJMYjQEJIiQKSjXYY2cft6VhnSMwsV8sAKC0qgZpx4sb9PWSJKHGjSfRE/kSZqQTSBLw3HPAN9+IWYu//up7Q+vdQEPzEXBMRgKAn1KB2ZOSEKDywy/7C5B9tqz+L4L35SMPoPAiM2bMwPnz5wEAY8aMwV133QWdTgc1h+Q1zpEj4njX4mJg0iRg6VLfmtzuRtYcOI9Xf8hGhdb+O1kBKiWGdY7E3D/1g1GScPPsrThbav+TJQBQ+ynQO7Y5vno4GWGBatw9d6fdoRMaoMJHE3vhtsToBn0mETlH7YwcM2YM/Pz8oFRy+3ijcNaiW2hMPgJNz8hAtRL/eqAfhneNwqajF/HnxVmoMdhXIHlbPvLJlBeZ8Uf7mUajQVBQEKqqqvDMM8/g/vvvx6233irv4jxNXh5wyy2ikLrtNmDxYhZSMgpUK6GwY35GgEoJncEIowRo9UZsyylCZY0eYYFqxLYMtisoVEoFlAoFagxG6AwSMk9dRmG5FmGBaiREhNhdTCkABKr4Z4bIXVhmpEKhQE5ODl5//XV8//33CORA2YaZOVMUUioV8P33LKRkZG8+Ao7JSEDsJf764WT069ACP2Sdxas/7IehAW2C3paPLKa8kCkUFi9ejIULF+J///sfdu3aha5du8q8Mg9RUACMHAmcOyeeTK1YAdh5khw5R3xEqF393H1im+Mf9/VBRl4JyjQ6JMe1RFigGsUVWhw8J9oPJvZth9mTknC4oAxjvthh8x5twgLxy3NDsPPkJVyqrEFiTBg6RYVCU2NAxqkSu9esN0qI95LNtUTeJDAwEFqtFuPGjcOJEycwdepULFq0yGuOaXa6OXPE0HqFQrT4cWi9rOzNR8AxGQkAXz3UH/06tEBpVQ3KNTq8OaYbAGBbThG25RTVuw5vy0cWU17IFAhTpkzBunXrsHLlStx+++1IT0+/uuGW6lBcLJ5I5eYC/fsDq1YBwcFyr8rndWgZDKMdHckXr2iRV1yJIZ0jEOKvQkllDX7ZX4DPNx/HFa2YHG/6eamu8Lmi1SP7bCmS41oiPEiNMk0NNh8txJytJxp0lKtRktChFf/sELkbhUKBwMBArFixAoMHD8bixYvRvXt3TJs2Te6lub9Fi8TQegCYN49D692AvfkIOCYjAaBNuLhp3zzYH48Ojr/66+UanV3FlLflI4spL6ZUKrFo0SKcOnUKWVlZuPvuu7F+/Xr48ynLtZWXA6NHA4cPAz16iM20PC5XfqWlUK5fj85VNTjo3+K6l+YVV1od+3otXduITedzt5285utlGh0e/k9m49ZqoUvrUN7pJnJjSUlJWLJkCcaPH4833ngDN9xwAyZMmCD3stwXZy26H0mC8tBBdDZW4iCC6r3cERkJAEM+3NKwddbibfnIXZdeLiQkBD/99BOio6Oxbds2PPvss+CZI9dQVSVaFfbsARISxLwMPsWThyQBBw8CH34I3HgjEBEB3HsvkrO3Q2Fs+uk/gzpG4Ofsc1hzsP5ZGI2lUADJcfzzQ+Tu7rrrLsyaNQsAMHnyZOzdu7eer/BRnLXoPiorgZ9/Bp56CujQAejVC8np6x2Sj4DzM9Ib85Gn+fmIzMxMDBs2DNXV1fj000/x4osvyr0k91FTA9x1l3gS1bYtkJYGxMfX/3XkOFVVwObNwOrVwJo1wJkz5tf8/IAhQ5B2yyQ8VR2PSr37f8sK9vfD/Mn9MaRThNxLIaJ6SJKERx99FP/973/Rtm1bZGZmIjraO04Zc4jffhPt7xqNmLX46accEeJqJ0+KfFy9Gti6VfzcYtK6NdLGPYKnWg1BpdH9/7t4Yz6ymPIhy5cvx7333gulUolVq1ZhzJgxci9Jfno9cP/94pCJiAhRSPGgDtfIzTUXT1u2AFqt+bWoKHGK4tixIsSbN4fRKGHArE0ovKKt+z3dROtmAfj99RFQNmAqPRHJR6vVYsSIEfjtt9+QnJyMbdu2ISio/rYpr5eVJU7qKy8Xw3kXLAB4lLzzabXi55E1a0RO5uSYX1MogJQUMQNz7FigTx8YoWA+yojFlI955513MGPGDDRr1gy///47evToIfeS5GM0AlOmAP/5DxAeLn6g79NH7lV5r5oa63A4dsz69eRkEQxjxgD9+l0zsBek5eLjDcdQrXPfYX9BaiVevuUGTBmaIPdSiKgBioqKkJKSglOnTuHee+/F0qVLvWpfR4Nx1qJrnTsn8nHNGmDjRqCiwvxa8+bAqFEiI0eNEjcca2E+yofFlI+RJAn3338/li1bhvj4eKSnpyMyMlLuZbmeJAEvvgh8/rk4rW/9emDwYLlX5X0KCsQwx9WrxT40y3AIDxehMGaMOPijdet6365Mo0PKBxuhdePJ6QEqJTLeGInwIA7LJvI0Bw8exMCBA1FRUYF33nkH06dPl3tJ8sjLA4YOFT/g33Yb8OOPHBHiaAYDkJ5ubt/LzrZ+vWdP8w3GgQPFTK/rYD7Kh8WUD9JoNLjxxhuRmZmJIUOGYOPGjQgICJB7Wa41fTrw3nsiHH75RbSSUdMZDEBGhrl9r/Zm7sREc2vCwIGAuuHfUKf/dBDLd+ej2g0DI1ClxD3JsXj3zkS5l0JEjbR69WrccccdkCQJy5Ytwz333CP3klyroEAUUrm54snU2rUcEeIoly6J/dlr1oi/l1jMLgwOFjMux4wRBWz79g1+e+ajPFhM+aiCggKkpKTg3LlzeOSRR/D111/7TjvD7NnAX/8q2hW+/x4YP17uFXm2khJg3TpRQP36qwgLk6AgYMQIUTzddps4eaiJNDUG3Dh7i1v2hrcOC8DWl29GkD9bYYg82SeffIKXX34ZgYGB2L59O5KTk+VekmsUF4tTVA8fFrMWN23iiJCmkCRg3z5ze3t6uthiYNKxo8jHsWNF4RoY2KSPYz7Kg8WUD8vKysLQoUNRVVWFDz/8EH/961/lXpLzzZ8PPPmk+OdFi4DJk+VdjyeSJGD/fvPTp99/tw6HhARza8JNNzU5HK4l81QJJn+d7la94YFqJRY/lor+cS3lXgoRNZEkSZg6dSq++uorREdHIzMzE23btpV7Wc5VXg4MHy5GhPToAWzbxhEhjXHlitjzZMrI8+fNr6nVolg1ZWSXLg7/eOaj67GY8nErV67E3XffDYVCgR9//BF33nmn3EtynqVLgQcfFMXAnDnA00/LvSLPUVEhwsG0OfbcOfNrarW4o2Zq3+vSxSXH5r714wGs2HPWLdoZAlVKTOwXi7+P8772BSJfVVNTg1tvvRXbtm1D3759sX37doSEhMi9LOeoqhLdA9u3i6claWkAj4e3jySJ0/ZMxdP27YBOZ349JsacjyNGAM2aOX1JzEfXYjFF+OCDD/Dmm28iJCQEO3fuRK9eveRekuOtWiXa+QwGYOZM4PXX5V6R+8vJMbcmbN9uPdciOto6HGRoA9HqDbh/wS4cKiiXdcNtgEqJHjFhWDp1AAJU3te+QOTLLl26hNTUVJw8eRJ33303li9fDqW3HQ1ee9bijh1AXJzcq3Jv1dVi3pMpI3Nzza8plWJPsCkje/Vy+Vwu5qNrsZgiSJKEyZMnY8mSJWjfvj0yMjLQ2o6T1TzG5s3im5pWK4qomTPlXpF7qq4WRZPp7tqJE+bXFApgwABza0Lv3m4xtLFSq8fEeTuRW1QpS2AEqJRIiAzBiicHISTg+ictEZFnOnLkCAYMGIDy8nK89dZbeO+99+RekuNw1qL9zpwxF0+bNokhxiatWolTaceOBW691S3aI5mPrsNiigAA1dXVuPnmm7Fr1y4MHDgQmzdvRqAT9rq43K5d4nScykrR1vfPf7pFEeA28vOt51pUVZlfa9lShMOYMeII8wj3nFZeqdVj8tfpOFxQ7tKWhkCVEt1jwvDNY6leHxREvm7dunUYM2YMjEYjlixZggceeEDuJTUdZy1en04n9gSbbjAePGj9ep8+5qdPKSluOYOL+egaLKboqosXLyI5ORn5+fn405/+hEWLFnn2CX/794uNnqWl4qCJhQs5uV2vF+Fgurt24ID16717i3AYMwZITa13roW70OoNeO+XI1iRle+STbeBaiUm9o3F27d38+rWBSIy++KLL/D8888jICAAW7duxYABA+ReUuNx1uK1FRaKo+DXrBGn1JaVmV8LDRVjVEyn08bEyLfOBmA+Oh+LKbKSnZ2NwYMHo7KyEu+//z7eeOMNuZfUODk5Yk5GYaHYK7V8uccUBg2h1WrrnxFWVCR64VevFuFQWmp+LSTEOhw8/LSqzFMleObbLJRrdE65CxeoUiIsSI05D/RFspeeSkRE1yZJEp5++mnMnTsXrVu3RkZGBto3YhaQW/CBWYt25aPRKE4vNN1g3L1bFJomN9xgbm8fOtSjBxczH52HxRTZ+PnnnzFu3DhIkoQffvgBEyZMkHtJDXPmDDBkiGhhu+UWcfiEFw0lvnjxIubOnYvly5fjnnvuwSuvvGJ9wpTRKIblmloTMjKsw6FLF+tw8KJ/N4CYszHz1yNYlpkPpQLQOOBOXJBaCaME3Jsci2mju3nlnAwiqp9Op8Po0aOxefNmJCUlYceOHQgNDZV7WQ3jxbMW681HQNxQ3LBBZOTateKmq0lAAHDzzeYOjY4dXbp+Z2M+OgeLKbqmDz/8EK+99hqCg4ORlpaGvn37yr0k+1y8KAqE48eBQYNE64KXHWW7aNEibNy4Ec8//zx69+4Ng8FgvvtmNAIrVwKTJpm/wN/fOhw6dZJn4S5WptHh+935mJ+WiwqtHhqdAQ35bqdQAEFqP4QGqPDE0ARM6h+L8CC18xZMRB6hpKQEAwYMwPHjx3HXXXdh5cqVnnPCn5fPWrxuPgLioKW+fYEjR8y/FhtrHpx7881e9zPDtTAfHYvFFF2TJEl47LHHsHDhQrRt2xaZmZmIdveZEyUlYkjsgQNi78+WLUDz5jIvyrEuX76MBx98EN999x3C6jqO/PJlcRSr6enTiBE+EQ51MRol/HayGFuOFSEj7xKOF1ZAqVBApVRAgvizrlAooACgN0owShK6tA5FclwrDO8ahUEJraBUevDeQSJyuJycHKSmpqK0tBSvvfYaZs2aJfeS6uflsxbtyseqKuDtt0U7nykje/Tw2YOpmI+OwWKK6qTVajFy5Ejs2LEDycnJ2LZtG4KCguRe1rVduSJa+tLTxbGu27cDkZFyr8opevXqhTfeeAMbNmzAhQsX8PTTT2PYsGFoZhoEaDCIgzZ8NBzqYzRKOFNShbziSlTrDajRG+GvUiJQ5Yf4iBB0aBXs2QevEJFLbNq0CaNGjYLBYMDChQvx8MMPy72kuvnIrMV68xEQBzF54R5qR2A+Ng6LKbquoqIipKSk4NSpU7j33nuxdOlS9/sfqbpa3F3asgXo0EEMHGzXTu5VOUZJiRiI+8c3fr1ej5deegnZ2dlXDwdZsmQJevXqhVdeeUXOlRIR+Zy5c+fiz3/+M9RqNTZv3owhQ4bIvSRblrMWp00DPvhA7hU5Rk2NuHHIfCSZsTSn64qMjMSqVaswaNAgLFu2DN27d8f06dMb/D5Go4TTJVXIK65Atc4IncEItZ8SgWol4iNC0aFlcOMeFet0wD33iEIqOloM0vPkQkqnA377zXx4hCQBmZlXw0KlUqF9+/ZYtWoVRo0aBaPRiOrqaixdulTmhRMR+Z6nnnoKhw8fxhdffIHx48cjMzMTcXFxDXoPp+UjIGYt3nmnKKSeeQZ4//3GvY+7OH9eHBqxerU4RGLePODeewGlkvlIsuGTKbLLmjVrcMcdd8BoNGLZsmW45557rnu9qQ9389FCZJ4qqdWHK0GSRBeaAoqrfbido0KRHNcSw7tGYXDHiPrDw2AQm2eXLhUDZrdvF73PnubCBfNci/XrgfJy82thYeJ0wvDwq79UWFiIwYMHY8+ePQgLC8N7772H8PBwPPvss56zCZqIyEvo9XqMHTsW69evR2JiIn777be69+zARfkIeMesRYNB3FA03WDMyrJ+/cUXRdtiYCAA5iPJg8UU2e3TTz/FSy+9hMDAQGzfvh3Jyck215RpdFi+Ox8L0nJRqdWjqsaAhvwBUwAI9vdDSIAKU4cm4J66ToiRJHEi0YIFQLNm4onUNdbjloxGEQ6muRZ79li/3r27ear64MGA2vb3/8EHH+DMmTPIzs6Gn58fPv/8c885cZGIyMuUlpZi4MCBOHr0KMaOHYuffvoJfn7WR0S7LB8Bz561WFIiZiKuWSNmJBYXm18LChKHKplOp+3QwebLmY/kaiymyG6SJOGJJ57Al19+iejoaGRkZKDdHy11mhoDZq49gmW786FQwCFTtq/OLugfi2m3WcwukCQxI+Pjj8XdqHXrgGHDmvx5TnX5snjqtHq1CIeiIvNrgYHA8OHmwbnx8Xa9ZXZ2NrRaLVJSUpy0aCIisteJEyeQmpqKkpISvPTSS/j4448BuDgfAc+btShJ4ima6Qbj77+Lm44m8fHmo8tvvFEUVPVgPpIrsZiiBqmpqcGtt96Kbdu2oW/fvti+fTsOFWrx7FIXTtV+7z0xvV2tBn76SRQg7kaSgIMHza0JO3eKdgWTuDjzsaw332xXOBARkXvbunUrbrnlFuj1eixYsAC9RkxwbT56yqzFigrRUWLKyHPnzK+pVOIGqSkjb7iBp9OSW2MxRQ126dIlpKam4uSp0+g/dSauRPZwyJ22+gSqlZiouoy3330YAZIR+O476+G0cqusFOGwZo34Kz/f/JpKJQLO1L7XtSvDgYjIC3311VeY8uRTaDXySbTsdxtqnB+PIh97ROLtmVMRkL3PPWctHj9uLp62bROn8ZlER5tb90aOFPuFiTwEiylqlL0HDuOOT9ZB2aoDlGrXtQ8E6rTofvEkvrkxAiFTHnXZ59bpxAlza8LWrdbh0Lq1uXgaOdLqEAkiIvJOlVo9hk5fjmJDoGvz0aBD9/PH8c3exQjZslH+WYtarSiaTBl54oT5NYUCGDDAnJG9e/MGI3ksFlPUYJVaPSbO24kTF65AJ8OfngAYkRAdjhVPDkJIgIs31Gq1QFqa+e5aTo75NYUCSEkxtyb06eN5JycREVGjmfIxt6gSWie09dUnQF+DhNZhWPHcja7PRwA4e9ZcPG3aJDo2TFq0AEaPFhk5ahQQEeH69RE5AYspahCt3oD75u/C4fPlsgSFSYBKiR4xYVg6dQACVH71f0FTnDtnbt3buFH0eps0by7CYcwY8Xe57wQSEZEsfDIf9Xoxy8p0g3H/fuvXk5LMNxhTUz3nREGiBmAxRQ3y1o8HsGLPWadspG2oQJUSE/vF4u/jEh37xgaDCAfT3bXsbOvXe/UytyYMGMBwICIi38hHQJxG++uvIiPXrROn1ZqEhIgTBE37n9q2dfznE7kZ/hRIdsvIK8GKLPcICgCo1huxIisfd/WOEacYNUVxsQiF1avF30tKzK8FB4s9T6ajy2Njm/ZZRETkVbw6H41GYN8+89On9HRxYq1Jly7mG4xDh7r3MexETsAnU2QXTY0BN87egsIrWrmXYiOqWQC2vXKz9ZyN+kiSCAfT06f0dOu5Fp06mVsThg27Ol2diIjIktflIwCUlwMbNphb3C9cML/m7w/cdJM5Izt1cuiaiTwNn0yRXT5YewTlGp3cy7imco0OM389gnfvrKed4coVsefJdHft/Hnza/7+Yhig6e5a587OXTQREXkFr8hHSQKOHjXfYExLE/uhTNq1MxdPI0a45+wqIpnwyRTVq0yjQ8oHG2XdUFufAJUSGW+MRHiQ2vyLkiRO2zMVT9u3AzqLwGvb1lw8jRgBhIa6fuFEROSxPDYfAUCjESM9TBmZl2d+zc9PDP01ZWRiIo8uJ6oDn0xRvZbvznf776FKBfD97nxMSY4R4WC6u5aba3GREhg82Hx3rVcvhgMRETWaR+Xj0ATg9Glz8bR5syioTCIixL7gsWOBW28VR5kTUb34ZIquy2iUMGDWJrfsBa8tSleJXf96FEpNlfkXW7US4TBmjJhr0bKJG3GJiIjgYfkoabFr1XQoDx+yfqFfP/MNxv79xRMpImoQPpmi6/rtZDEqtfr6L3QDFZISO6M6Y0grP3NrQnIyw4GIiBzOo/JRZ8TOShWGhIWJp05jxogbjW3ayL00Io/HYoqua/PRQlTVGORehl00/oHYPGs+htyXIvdSiIjIy3lUPqoDsfn1/8OQx4cDanX9X0BEdmMxRdeVeaoE9vSBzhzfE/07tEBM8yDUGIzYl1+KD9YcwfHCCrs/S6EAXhjeGfcmx6JliD9OFlbgw/XHsPVYkV1fL0GBzGL3b7cgIiLPZ28+Ao7JyFdH3YA7k2IQGRqAar0Rxy5cwacbc/B77qV6v1ZSKJCpD2YhReQESrkXQO7LaJTs/kZ/f0p7VGj1+Dm7ABXVetx8QxQWPZaCAJX9f8SeGtYRL47sAr1Bwi/7z6NjZCi+nNwfnaPsP2Uv52IFuA2QiIicqSH5CDgmI2NbBCP7bBmW7zmLMyVVSIlvia8fTkaQ2r5WduYjkXPwyRTV6XRJFZR2HlM04d87kXXmMgCgXfMg7HhtOKLDg9ApKhSHCsrr/Xo/pQJThyYAAP68ZA8OFpTjXKkGzw/vjCeHJeCVFfvtWodSocDpS1WIi+AMDCIico6G5CPQ9IwEgOe+23v1n8OD1MiefiuC/P0QEeqP/Mua63ylwHwkcg4WU1SnvOIKqJT2hYUpJABA/cedNoNRsvuUo+jwQLQM8YfBKOHgH8Fy4GwZAKB7dJjda1YpFcgrrmRYEBGR0zQkH4GmZ6TJnUkx6NehBfq2F8eW/7K/wK5CCmA+EjkLiymqU7XOCMnujnAh2N8PsycmAQAWpOWiyM6giAwNAABodObNvFU14pSkyGYBdn++BKBa7xkbgomIyDM1Jh+BxmekybDOEZjYLxYAUFpVg7TjxXZ/LfORyDm4Z4rqpDMY0ZD26hbBanw7ZQD6dWiBbzPOYNavR+3+2qIKEShBar+rAxBDAkSt35CwkSQJNW48iZ6IiDxfQ/MRaFpGmryyYj86v7UGk+bthFKpwP/d3Qv9Otg3XJf5SOQcLKaoTmo/pd2T3ds2D8KKpwahd2xz/GvrCbzxvwMN+qzzZdW4XFUDP6UCPduGAwB6tWsOADhy4Yrd76NQKODfgA29REREDdWQfASanpEqpQL+fiLbdAYJmacuo7Bc3GhMsLNtj/lI5Bxs86M6BaqVUMC+tPjhqUFoEx6Is5erEKj2w/TbuwMAftp3DtlnyzCxbzvMnpSEwwVlGPPFDpuvNxglLEjLxaujuuJfD/RFel4Jbu8ZDb3BiHnbT9q9ZgWAQBWH9BIRkfM0JB+Bpmdkm7BA/PLcEOw8eQmXKmuQGBOGTlGh0NQYkHGqxK41MB+JnIPFFNUpPiIUeqN9fQxtwgMBAO1aBOOxwfFXf/1wQTmyz5ZdvYN3vfebu+0kAtV+uKdfLO7oFYOTRRX4aP0x5Fy0//hZvVFCPDfXEhGREzUkH4GmZ+QVrR7ZZ0uRHNcS4UFqlGlqsPloIeZsPYHTl6rsWgPzkcg5WExRnTq0DIbRzqbwuGmrr/t61zbNAIiCqS5GCfhkQw4+2ZBj/yJt3kNCh1bBjf56IiKi+jQkH4GmZ2SZRoeH/5Np/wKvgflI5BxsnqU6KZWKBg3MvZ5BHSPwc/Y5rDl4wSHvV5curUOhaEgjOxERUQM5Mh8B12Qk85HIOfhkiq4rOa4lDhWUN+IAWGu3fZ7mkPVcj0IBJMe1cvrnEBEROSofAednJPORyHn4ZIqua3jXKAT7e8aG1SC1H4Z3jZJ7GURE5AOYj0QEsJiiegzuGHF13pO7axagwqAE3nkjIiLnYz4SEcBiiuqhVCowdWgCAtXu/UclSK3E1KEJUCrZD05ERM7HfCQigMUU2eGe/rENnvTuakYJmNQ/Vu5lEBGRD2E+EhGLKapXeJAa9/aPRaCbTk4PVClxb3IswoPUci+FiIh8CPORiNzz/35yO9Nu64YwN/1mHB6sxrTR3eReBhER+SDmI5FvYzFFdgny98OcB/q6XW94oFqJOff3RZCHnKhERETehflI5Nvc6/98cmvJcS0xsW87t2lnCFQpMbFvLPrHtZR7KURE5MOYj0S+yz3+ryeP8fbt3dE9JgwBMgdGgEqJ7jFhePt2ti8QEZH8mI9EvonFFDVIgMoP3zyWioTIENkCI0ClREJkCL55LBUBKrYvEBGR/JiPRL6JxRQ1WEiACiueHIQeMWEub2kIVCnRIyYMK54c5DHDEomIyDcwH4l8j0KS3H1CArkrrd6A9345ghVZ+ajWGZ3+eYFq0QP+9u3deMeNiIjcFvORyHewmKImyzxVgme+zUK5RodqveNDI1ClRFiQGnMe6ItkbqYlIiIPwXwk8n4spsghNDUGzPz1CJZl5kOpADQOuBMXpFbCKAH3Jsdi2uhuPN6ViIg8DvORyLuxmCKHKtPo8P3ufMxPy0WFVg+NzoCG/AlTKIAgtR9CA1R4YmgCJvXn5HYiIvJ8zEci78RiipzCaJTw28libDlWhIy8SzheWAGlQgGVUgEJgCRJUCgUUADQGyUYJQldWociOa4VhneNwqCEVlAqFXL/NoiIiByK+UjkXVhMkUsYjRLOlFQhr7gS1XoDavRG+KuUCFT5IT4iBB1aBUOhYDgQEZFvYT4SeTYWU0RERERERI3AOVNERERERESNwGKKiIiIiIioEVhMERERERERNQKLKSIiIiIiokZgMUVERERERNQILKaIiIiIiIgagcUUERERERFRI7CYIiIiIiIiagQWU0RERERERI3AYoqIiIiIiKgRWEwRERERERE1AospIiIiIiKiRmAxRURERERE1AgspoiIiIiIiBqBxRQREREREVEjsJgiIiIiIiJqBBZTREREREREjcBiioiIiIiIqBFYTBERERERETUCiykiIiIiIqJG+H8vByUxbIv5IAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1080x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "label_dict = {i: str(i) + \", \" + str(t) for i, t in solution.items()}\n",
    "edge_color = [\"red\" if solution[u] == (solution[v] + 1) % n\n",
    "              or solution[v] == (solution[u] + 1) % n else \"black\"\n",
    "              for (u, v) in G.edges]\n",
    "label_dict_bf = {i: str(i) + \", \" + str(t) for i, t in solution_brute_force.items()}\n",
    "edge_color_bf = [\"red\" if solution_brute_force[u] == (solution_brute_force[v] + 1) % n\n",
    "                 or solution_brute_force[v] == (solution_brute_force[u] + 1) % n else \"black\"\n",
    "                 for (u, v) in G.edges]\n",
    "\n",
    "# 在图上画出上面得到的最优路线\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 4))\n",
    "for i, a in enumerate(ax):\n",
    "    a.axis('off')\n",
    "    a.margins(0.20)\n",
    "nx.draw(G, pos=pos, labels=label_dict, edge_color=edge_color, ax=ax[0], **options)\n",
    "nx.drawing.nx_pylab.draw_networkx_edge_labels(G, pos=pos, ax=ax[0], edge_labels=nx.get_edge_attributes(G, 'weight'))\n",
    "nx.draw(G, pos=pos, labels=label_dict_bf, edge_color=edge_color_bf, ax=ax[1], **options)\n",
    "nx.drawing.nx_pylab.draw_networkx_edge_labels(G, pos=pos, ax=ax[1], edge_labels=nx.get_edge_attributes(G, 'weight'))\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面给出的左图展示了参数化量子电路找到的旅行商问题的最优解，右图展示了经典暴力算法找到的最优解。我们不难看出，即使旅行商访问每个城市的绝对顺序不一样，但路线是一致的，即相对顺序一样。这说明在这个例子中，参数化量子电路找到了旅行商问题的最优解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  实际应用\n",
    "\n",
    "旅行商问题可以直接应用在很多交通和物流规划中，例如规划校车接送学生的路线。管理科学领域的先锋 Merrill Flood 在上世纪40年代就受到校车问题的启发，展开了对于旅行商问题的早期研究。更多近期的应用包括了路线规划 [1] 和电线公司对于电力传输的规划 [2]。\n",
    "\n",
    "除了交通运输问题，旅行商问题同样在管理问题中有很广泛的应用，比如计划机器在电路板上钻孔的顺序 [3]、重构 DNA 上的不明片段 [4] 以及规划最佳建筑路线 [5]。 一些咨询公司，比如 [Nexus](https://nexustech.com.ph/company/newsletter/article/Finding-the-shortest-path-Optimizing-food-trips) 也通过旅行商问题给外界提供管理咨询服务。\n",
    "\n",
    "同时作为最著名的组合优化问题之一，旅行商问题为很多用于解决组合问题的通用算法提供了测试平台。它经常被作为研究者测试他们提出的新的算法的首选例子。\n",
    "\n",
    "对于旅行商问题更多的应用和解法，详见 [6]。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Bräysy, Olli, et al. \"An optimization approach for communal home meal delivery service: A case study.\" [Journal of Computational and Applied Mathematics 232.1 (2009): 46-53.](https://www.sciencedirect.com/science/article/pii/S0377042708005438)\n",
    "\n",
    "[2] Sloane, Thomas H., Frank Mann, and H. Kaveh. \"Powering the last mile: an alternative to powering FITL.\" [Proceedings of Power and Energy Systems in Converging Markets. IEEE, 1997.](https://ieeexplore.ieee.org/document/646046)\n",
    "\n",
    "[3] Onwubolu, Godfrey C. \"Optimizing CNC drilling machine operations: traveling salesman problem-differential evolution approach.\" [New optimization techniques in engineering. Springer, Berlin, Heidelberg, 2004. 537-565.](https://link.springer.com/chapter/10.1007/978-3-540-39930-8_22)\n",
    "\n",
    "[4] Caserta, Marco, and Stefan Voß. \"A hybrid algorithm for the DNA sequencing problem.\" [Discrete Applied Mathematics 163 (2014): 87-99.](https://www.sciencedirect.com/science/article/pii/S0166218X12003253)\n",
    "\n",
    "[5] Klanšek, Uroš. \"Using the TSP solution for optimal route scheduling in construction management.\" [Organization, technology & management in construction: an international journal 3.1 (2011): 243-249.](https://www.semanticscholar.org/paper/Using-the-TSP-Solution-for-Optimal-Route-Scheduling-Klansek/3d809f185c03a8e776ac07473c76e9d77654c389)\n",
    "\n",
    "[6] Matai, Rajesh, Surya Prakash Singh, and Murari Lal Mittal. \"Traveling salesman problem: an overview of applications, formulations, and solution approaches.\" [Traveling salesman problem, theory and applications 1 (2010).](https://www.sciencedirect.com/topics/computer-science/traveling-salesman-problem)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.13 ('new_pq')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "58b83104798ee1b81625bc6249d4d66f2bacd4a7d411a9b7a27bac0b4765adf2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
